In this section, we import the hydrological data and perform initial processing, which includes checking for any missing values and cleaning the dataset. If missing data is found, appropriate methods will be used to address them.
# Import dataset from excel file
main_data = read_excel("Savinja-VelikoSirje.xlsx")
# Create dataframe
main_data = as.data.frame(main_data)
# Convert the 'Date' column to Date type
main_data$Date <- as.Date(main_data$Date, format = "%d.%m.%Y")
head(main_data)
## Date Year Month Day Waterlevel Discharge Watertemp P1_Luce P2_Solcava
## 1 2010-01-01 2010 1 1 251 59.058 6.9 7.0 4.8
## 2 2010-01-02 2010 1 2 246 53.925 7.1 3.0 0.6
## 3 2010-01-03 2010 1 3 239 47.681 6.6 0.1 2.7
## 4 2010-01-04 2010 1 4 231 41.040 5.3 0.0 0.0
## 5 2010-01-05 2010 1 5 228 38.694 4.4 3.3 2.0
## 6 2010-01-06 2010 1 6 226 36.966 4.2 6.2 4.6
## P3_Lasko P4_GornjiGrad P5_Radegunda Airtemp Evapotrans
## 1 4.3 4.6 2.2 5.7 0.3
## 2 2.1 3.2 2.0 3.5 0.9
## 3 1.5 0.0 0.1 -2.0 0.5
## 4 0.0 0.0 0.0 -3.0 0.2
## 5 7.0 3.1 4.0 -1.8 0.2
## 6 7.0 9.2 8.2 -1.0 0.2
# Structure of dataframe
str(main_data)
## 'data.frame': 4748 obs. of 14 variables:
## $ Date : Date, format: "2010-01-01" "2010-01-02" ...
## $ Year : num 2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
## $ Month : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : num 1 2 3 4 5 6 7 8 9 10 ...
## $ Waterlevel : num 251 246 239 231 228 226 223 220 223 223 ...
## $ Discharge : num 59.1 53.9 47.7 41 38.7 ...
## $ Watertemp : num 6.9 7.1 6.6 5.3 4.4 4.2 4.2 4.2 3.7 3.6 ...
## $ P1_Luce : num 7 3 0.1 0 3.3 6.2 2.5 2 37.6 2.4 ...
## $ P2_Solcava : num 4.8 0.6 2.7 0 2 4.6 0.4 2 32.2 4.5 ...
## $ P3_Lasko : num 4.3 2.1 1.5 0 7 7 8.1 3.1 20.1 3.2 ...
## $ P4_GornjiGrad: num 4.6 3.2 0 0 3.1 9.2 3.9 2.8 36.1 4.3 ...
## $ P5_Radegunda : num 2.2 2 0.1 0 4 8.2 2.8 1.5 34.1 3.8 ...
## $ Airtemp : num 5.7 3.5 -2 -3 -1.8 -1 -0.3 -0.4 0.6 0.9 ...
## $ Evapotrans : num 0.3 0.9 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
# Check for missing values
missing_data <- colSums(is.na(main_data))
missing_data
## Date Year Month Day Waterlevel
## 0 0 0 0 183
## Discharge Watertemp P1_Luce P2_Solcava P3_Lasko
## 0 224 0 0 184
## P4_GornjiGrad P5_Radegunda Airtemp Evapotrans
## 0 0 0 0
sum(is.na(main_data$Waterlevel))
## [1] 183
# Fit linear regression model using Discharge
model <- lm(Waterlevel ~ Discharge, data = main_data, na.action = na.exclude)
model
##
## Call:
## lm(formula = Waterlevel ~ Discharge, data = main_data, na.action = na.exclude)
##
## Coefficients:
## (Intercept) Discharge
## 173.35 0.88
# Predict missing Waterlevel
missing_waterlevel <- which(is.na(main_data$Waterlevel))
predicted_values <- predict(model, newdata = main_data[missing_waterlevel, ])
predicted_values
## 2009 2010 2011 2012 2013 2014 2015 2016
## 197.4418 195.2383 193.6437 192.6052 191.6953 190.5953 189.9546 217.9093
## 2017 2018 2019 2020 2021 2022 2023 2024
## 203.1065 195.1661 192.8129 191.1277 190.0981 188.8528 188.0775 187.6472
## 2025 2026 2027 2028 2029 2030 2031 2032
## 186.7047 186.0421 185.6117 185.0608 185.3196 184.5997 188.0907 201.2136
## 2033 2034 2035 2036 2037 2038 2039 2040
## 197.9663 196.5495 200.5580 194.4436 227.3713 227.5033 203.9716 197.2421
## 2041 2042 2043 2044 2045 2046 2047 2048
## 196.7756 195.4697 192.1256 190.0514 188.7631 187.7326 186.7461 186.1406
## 2049 2050 2051 2052 2053 2054 2055 2056
## 185.2861 184.8743 184.1685 183.8543 183.8367 184.4607 186.8376 186.7883
## 2057 2058 2059 2060 2061 2062 2063 2064
## 187.9411 205.1640 187.8223 185.2923 184.1333 183.6141 183.5384 187.5610
## 2065 2066 2067 2068 2069 2070 2071 2072
## 186.4134 184.5117 183.6141 183.0482 182.6170 182.1700 181.9465 182.0213
## 2073 2074 2075 2076 2077 2078 2079 2080
## 184.0268 222.2223 237.3929 208.2124 196.2529 191.6548 189.4240 187.9271
## 2081 2082 2083 2084 2085 2086 2087 2088
## 186.8904 186.0421 185.4190 198.7214 193.9270 190.3366 188.5862 187.6692
## 2089 2090 2091 2092 2093 2094 2095 2096
## 186.9133 185.9030 185.3433 184.9491 193.3603 218.6925 220.0477 202.5609
## 2097 2098 2099 2100 2101 2102 2103 2104
## 195.5586 191.9012 190.0285 188.8308 187.8611 187.3700 186.7866 187.7150
## 2105 2106 2107 2108 2109 2110 2111 2112
## 186.8904 229.9902 225.4678 207.6984 198.9590 264.2088 261.3857 228.6218
## 2113 2114 2115 2116 2117 2118 2119 2120
## 307.3579 510.5842 405.4642 279.2572 245.3932 254.7803 258.4201 232.4146
## 2121 2122 2123 2124 2125 2126 2127 2128
## 220.7342 213.7688 208.5626 205.1464 202.5169 200.2508 198.5858 197.0643
## 2129 2130 2131 2132 2133 2134 2135 2136
## 195.4847 194.1233 192.7856 191.8203 191.2755 190.7167 190.0972 189.4900
## 2137 2138 2139 2140 2141 2142 2143 2144
## 189.0104 188.4067 188.0775 187.7106 187.2635 186.8490 186.4812 186.2198
## 2145 2146 2147 2148 2149 2150 2151 2152
## 185.8045 185.8045 185.8837 185.1559 184.9508 185.6901 218.8606 224.1381
## 2153 2154 2155 2156 2157 2158 2159 2160
## 207.3922 198.0130 195.1529 193.7554 191.7763 190.3814 189.4451 188.7842
## 2161 2162 2163 2164 2165 2166 2167 2168
## 188.8291 188.8520 189.2603 188.8317 187.5407 186.9106 186.4011 186.2401
## 2169 2170 2171 2172 2173 2174 2175 2176
## 185.8247 185.6311 184.9667 184.8558 184.6728 184.2935 184.0268 183.8367
## 2177 2178 2179 2180 2181 2182 2183 2184
## 183.6123 183.6123 183.6123 183.6123 183.5454 183.5789 183.2603 183.0808
## 2185 2186 2187 2188 2189 2190 2191
## 182.8370 182.8212 182.8212 182.8212 182.8212 182.8212 182.5871
# Fill missing Waterlevel values
main_data$Waterlevel[missing_waterlevel] <- predicted_values
# Check Waterlevel again
sum(is.na(main_data$Waterlevel))
## [1] 0
missing_watertemps = which(is.na(main_data$Watertemp))
missing_watertemps
## [1] 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992
## [16] 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
## [31] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022
## [46] 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037
## [61] 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052
## [76] 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067
## [91] 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082
## [106] 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097
## [121] 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112
## [136] 2113 2114 2115 2116 2117 2118 2119 2120 2121 2122 2123 2124 2125 2126 2127
## [151] 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 2140 2141 2142
## [166] 2143 2144 2145 2146 2147 2148 2149 2150 2151 2152 2153 2154 2155 2156 2157
## [181] 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168 2169 2170 2171 2172
## [196] 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187
## [211] 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 2199 2200 2201
# Loop missing watertemp and fill
for (i in missing_watertemps){
current_year <- main_data$Year[i]
date_without_year <- format(main_data$Date[i], "%m-%d")
# Find Watertemp and Airtemp for same date in previous years
previous_data <- main_data %>%
filter(format(Date, "%m-%d") == date_without_year & Year < current_year) %>%
dplyr::select(Watertemp, Airtemp)
# Calculate mean of Watertemp and Airtemp
if(nrow(previous_data) > 0){
watertemp_mean <- mean(previous_data$Watertemp)
airtemp_mean <- mean(previous_data$Airtemp)
main_data$Watertemp[i] <- watertemp_mean + (main_data$Airtemp[i] - airtemp_mean)
}
}
# Check Watertemp
sum(is.na(main_data$Watertemp))
## [1] 0
missing_lasko <- which(is.na(main_data$P3_Lasko))
# Fill missing values with average of other stations
for (i in missing_lasko){
other_prec <- main_data[i, c("P1_Luce", "P2_Solcava", "P4_GornjiGrad", "P5_Radegunda")]
main_data$P3_Lasko[i] <- mean(as.numeric(other_prec))
}
# Check P3_Lasko
sum(is.na(main_data$P3_Lasko))
## [1] 0
# Check main_data for any missing values
colSums(is.na(main_data))
## Date Year Month Day Waterlevel
## 0 0 0 0 0
## Discharge Watertemp P1_Luce P2_Solcava P3_Lasko
## 0 0 0 0 0
## P4_GornjiGrad P5_Radegunda Airtemp Evapotrans
## 0 0 0 0
In this section, we perform basic hydrological analysis. This includes calculating descriptive statistics for key variables such as discharge, precipitation, evapotranspiration, and air temperature. We will also plot graphs for these variables over the entire period and for the year 2022. Additionally, we will calculate correlations between precipitation stations, water balance, runoff coefficient, and the ratio between evapotranspiration and precipitation.
# Calculate descriptive statistics for Discharge
discharge_stats <- data.frame(
Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
Values = c(
mean(main_data$Discharge),
sd(main_data$Discharge),
min(main_data$Discharge),
max(main_data$Discharge),
median(main_data$Discharge))
)
discharge_stats
## Statistic Values
## 1 Mean 38.48014
## 2 Standard Deviation 51.68336
## 3 Minimum 5.51500
## 4 Maximum 859.74400
## 5 Median 22.03150
# Summary statistics for Discharge
summary(main_data$Discharge)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.515 13.966 22.032 38.480 40.252 859.744
# Calculate descriptive statistics for Precipitation
precipitation_stats <- data.frame(
Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
Values = c(
mean(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
sd(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
min(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
max(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda)),
median(c(main_data$P1_Luce, main_data$P2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda))
)
)
precipitation_stats
## Statistic Values
## 1 Mean 3.925959
## 2 Standard Deviation 9.844119
## 3 Minimum 0.000000
## 4 Maximum 135.400000
## 5 Median 0.000000
# Summary statistics for Precipitation
summary(c(main_data$P1_Luce, main_data$p2_Solcava, main_data$P3_Lasko, main_data$P4_GornjiGrad, main_data$P5_Radegunda))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.00 0.00 3.85 2.10 135.40
# P1_Luce
mean(main_data$P1_Luce)
## [1] 4.235573
sd(main_data$P1_Luce)
## [1] 10.816
min(main_data$P1_Luce)
## [1] 0
max(main_data$P1_Luce)
## [1] 116.1
median(main_data$P1_Luce)
## [1] 0
# P2_Solcava
mean(main_data$P2_Solcava)
## [1] 4.23155
sd(main_data$P2_Solcava)
## [1] 10.60273
min(main_data$P2_Solcava)
## [1] 0
max(main_data$P2_Solcava)
## [1] 109.5
median(main_data$P2_Solcava)
## [1] 0
# P3_Lasko
mean(main_data$P3_Lasko)
## [1] 3.255281
sd(main_data$P3_Lasko)
## [1] 8.140033
min(main_data$P3_Lasko)
## [1] 0
max(main_data$P3_Lasko)
## [1] 106.7
median(main_data$P3_Lasko)
## [1] 0
# P4_GornjiGrad
mean(main_data$P4_GornjiGrad)
## [1] 4.119756
sd(main_data$P4_GornjiGrad)
## [1] 10.10802
min(main_data$P4_GornjiGrad)
## [1] 0
max(main_data$P4_GornjiGrad)
## [1] 135.4
median(main_data$P4_GornjiGrad)
## [1] 0
# P5_Radegunda
mean(main_data$P5_Radegunda)
## [1] 3.787637
sd(main_data$P5_Radegunda)
## [1] 9.278106
min(main_data$P5_Radegunda)
## [1] 0
max(main_data$P5_Radegunda)
## [1] 110.5
median(main_data$P5_Radegunda)
## [1] 0
# Calculate descriptive statistics for Air Temperature
airtemp_stats <- data.frame(
Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
Values = c(
mean(main_data$Airtemp),
sd(main_data$Airtemp),
min(main_data$Airtemp),
max(main_data$Airtemp),
median(main_data$Airtemp)
)
)
airtemp_stats
## Statistic Values
## 1 Mean 10.878265
## 2 Standard Deviation 8.140201
## 3 Minimum -13.800000
## 4 Maximum 28.900000
## 5 Median 11.300000
# Summary statistics for Air Temperature
summary(main_data$Airtemp)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -13.80 4.40 11.30 10.88 17.70 28.90
# Calculate descriptive statistics for Evapotranspiration
evapotrans_stats <- data.frame(
Statistic = c("Mean", "Standard Deviation", "Minimum", "Maximum", "Median"),
Values = c(
mean(main_data$Evapotrans),
sd(main_data$Evapotrans),
min(main_data$Evapotrans),
max(main_data$Evapotrans),
median(main_data$Evapotrans)
)
)
evapotrans_stats
## Statistic Values
## 1 Mean 2.264006
## 2 Standard Deviation 1.713055
## 3 Minimum 0.000000
## 4 Maximum 7.700000
## 5 Median 1.800000
# Summary statistics for Evapotranspiration
summary(main_data$Evapotrans)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.700 1.800 2.264 3.600 7.700
# Plot Discharge over the whole period
plot(main_data$Date, main_data$Discharge, type = "l",
xlab = "Date", ylab = "Discharge (m^3/s)", main = "Discharge Over Whole Period")
# Plot Discharge for year 2022
data_2022 <- main_data[format(main_data$Date, "%Y") == "2022",]
plot(data_2022$Date, data_2022$Discharge, type = "l",
xlab = "Date", ylab = "Discharge (m^3/s)", main = "Discharge For Year 2022")
# Plot Precipitation over whole period for all stations
plot(main_data$Date, main_data$P1_Luce, type = "l", col = "blue",
xlab = "Date", ylab = "Precipitation (mm)", main = "Precipitation Over Whole Period")
# All lines for other stations
lines(main_data$Date, main_data$P2_Solcava, col = "red")
lines(main_data$Date, main_data$P3_Lasko, col = "green")
lines(main_data$Date, main_data$P4_GornjiGrad, col = "purple")
lines(main_data$Date, main_data$P5_Radegunda, col = "orange")
# Add legend
legend("topright", legend = c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda"),
col = c("blue", "red", "green", "purple", "orange"), lty = 1)
# Plot Precipitation for year 2022
plot(data_2022$Date, data_2022$P1_Luce, type = "l", col = "blue",
xlab = "Date", ylab = "Precipitation (mm)", main = "Precipitation for Year 2022")
# Add lines for other stations
lines(data_2022$Date, data_2022$p2_Solcava, col = "red")
lines(data_2022$Date, data_2022$P3_Lasko, col = "green")
lines(data_2022$Date, data_2022$P4_GornjiGrad, col = "purple")
lines(data_2022$Date, data_2022$P5_Radegunda, col = "orange")
# Add a legend
legend("topright", legend = c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda"),
col = c("blue", "red", "green", "purple", "orange"), lty = 1)
# Calculate the average precipitation across all stations
main_data$Avg_Precipitation <- rowMeans(main_data[, c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")])
# Plot Average Precipitation over the whole period
plot(main_data$Date, main_data$Avg_Precipitation, type = "l",
xlab = "Date", ylab = "Average Precipitation (mm)", main = "Average Precipitation Over Whole Period")
# Plot Average Precipitation for year 2022
data_2022 <- main_data[format(main_data$Date, "%Y") == "2022",]
plot(data_2022$Date, data_2022$Avg_Precipitation, type = "l",
xlab = "Date", ylab = "Average Precipitation (mm)", main = "Average Precipitation For Year 2022")
# Plot Evapotranspiration over whole period
plot(main_data$Date, main_data$Evapotrans, type = "l", col = "blue",
xlab = "Date", ylab = "Evapotranspiration (mm)", main = "Evapotranspiration Over Whole Period")
# Plot Evapotranspiration for year 2022
data_2022 <- main_data[format(main_data$Date, "%Y") == "2022",]
plot(data_2022$Date, data_2022$Evapotrans, type = "l", col = "blue",
xlab = "Date", ylab = "Evapotranspiration (mm)", main = "Evapotranspiration For Year 2022")
# Plot Air Temperature over the whole period
plot(main_data$Date, main_data$Airtemp, type = "l", col = "red",
xlab = "Date", ylab = "Air Temperature (°C)", main = "Air Temperature Over Whole Period")
# Plot Air Temperature for year 2022
plot(data_2022$Date, data_2022$Airtemp, type = "l", col = "red",
xlab = "Date", ylab = "Air Temperature (°C)", main = "Air Temperature For Year 2022")
# Precipitation columns
precipitation_data <- main_data[, c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")]
# Calculate correlation matrix
cor_matrix <- cor(precipitation_data)
cor_matrix
## P1_Luce P2_Solcava P3_Lasko P4_GornjiGrad P5_Radegunda
## P1_Luce 1.0000000 0.8974698 0.7280284 0.9175150 0.8899014
## P2_Solcava 0.8974698 1.0000000 0.7226416 0.8618554 0.8479587
## P3_Lasko 0.7280284 0.7226416 1.0000000 0.7861900 0.8029485
## P4_GornjiGrad 0.9175150 0.8618554 0.7861900 1.0000000 0.9092337
## P5_Radegunda 0.8899014 0.8479587 0.8029485 0.9092337 1.0000000
# Plot correlogram
corrplot(cor_matrix, method = "shade", type = "full",
tl.col = "orange", tl.srt = 45,
title = "Correlogram of Precipitation Stations",
mar = c(0, 0, 2, 0)
)
The
correlation analysis for selected precipitation stations in the Savinja
River basin reveals high inter-station correlations, with values ranging
from 0.72 to 0.92. The strongest correlation is observed between P1_Luce
and P4_GornjiGrad (0.92), indicating similar precipitation patterns,
likely influenced by geographic or climatic factors. The weakest
correlation (0.72) is between P2_Solcava and P3_Lasko, which may reflect
localized variations in precipitation. The correlogram visually
highlights these relationships, with darker shades representing stronger
correlations. Stations such as P1_Luce, P4_GornjiGrad, and P5_Radegunda
exhibit particularly high inter-correlations, suggesting potential
redundancy in data collection across these locations.
# Catchment area in m²
catchment_area <- 1847000000
# Calculate Runoff
main_data$Runoff <- (main_data$Discharge * 3600 * 24 * 1000) / catchment_area
main_data$Discharge_mm <- main_data$Discharge * 3600 * 24 * 1000
# Calculate the water balance
main_data$Water_Balance <- main_data$Avg_Precipitation - main_data$Evapotrans - main_data$Runoff
# Calculated Runoff
main_data$Calculated_Runoff <- main_data$Avg_Precipitation - main_data$Evapotrans
# Plot measured and calculated runoff
yearly_runoff <- main_data %>%
group_by(Year) %>%
summarise(
Total_Measured = sum(Runoff),
Total_Calculated = sum(Calculated_Runoff))
plot(yearly_runoff$Year, yearly_runoff$Total_Measured,
xlab = "Date", ylab ="Runoff", type = "l", col = "red", main = "Measured and Calculated Runoff")
lines(yearly_runoff$Year, yearly_runoff$Total_Calculated, lty = 2, col = "green")
legend("topright", legend = c("Measured", "Calculated"), col = c("red", "green"), lty = c(1, 2))
The comparison between measured and calculated runoff for the Savinja River basin shows a general agreement in the annual trends. While the measured runoff (red line) closely aligns with the calculated runoff (green dashed line), slight deviations are observed in certain years, such as 2012 and 2020, which may be attributed to discrepancies in precipitation or evapotranspiration data. Overall, the calculated runoff provides a reasonable approximation of the measured values, validating the water balance approach.
# Group data by Year and calculate total runoff and total precipitation
annual_runoff <- main_data %>%
group_by(Year) %>%
summarise(
Total_Runoff = sum(Runoff),
Total_Precipitation = sum(Avg_Precipitation)
)
# Calculate the runoff coefficient for each year
annual_runoff$Runoff_Coefficient <- annual_runoff$Total_Runoff / annual_runoff$Total_Precipitation
annual_runoff
## # A tibble: 13 × 4
## Year Total_Runoff Total_Precipitation Runoff_Coefficient
## <dbl> <dbl> <dbl> <dbl>
## 1 2010 816. 1603. 0.509
## 2 2011 369. 1168. 0.316
## 3 2012 625. 1514. 0.413
## 4 2013 793. 1484. 0.535
## 5 2014 1037. 1892. 0.548
## 6 2015 546. 1205. 0.454
## 7 2016 685. 1435. 0.478
## 8 2017 632. 1514. 0.417
## 9 2018 687. 1314. 0.523
## 10 2019 665. 1512. 0.440
## 11 2020 551. 1445. 0.381
## 12 2021 664. 1403. 0.474
## 13 2022 476. 1151. 0.413
# Plot the runoff coefficient over the years
plot(as.numeric(annual_runoff$Year), annual_runoff$Runoff_Coefficient, type = "b", col = "blue",
xlab = "Year", ylab = "Runoff Coefficient", main = "Runoff Coefficient Values")
The runoff coefficient values over the years show fluctuations between 0.31 and 0.55. Higher runoff coefficients, such as in 2013 and 2014, indicate a larger proportion of precipitation contributing to runoff, potentially due to higher precipitation intensity or reduced infiltration. Conversely, lower coefficients in years like 2011 and 2020 suggest greater infiltration or retention of precipitation within the catchment. This variability highlights the influence of climatic and catchment conditions on runoff generation.
# Group data by Year and calculate total evapotranspiration and total precipitation
evapotrans_prec <- main_data %>%
group_by(Year) %>%
summarise(
Avg_Evapotrans = sum(Evapotrans),
Avg_Precipitation = sum(Avg_Precipitation))
# Calculate evapotranspiration to precipitation ratio
evapotrans_prec$Ratio <- evapotrans_prec$Avg_Evapotrans / evapotrans_prec$Avg_Precipitation
evapotrans_prec
## # A tibble: 13 × 4
## Year Avg_Evapotrans Avg_Precipitation Ratio
## <dbl> <dbl> <dbl> <dbl>
## 1 2010 789. 1603. 0.492
## 2 2011 862. 1168. 0.739
## 3 2012 887. 1514. 0.586
## 4 2013 819. 1484. 0.552
## 5 2014 784. 1892. 0.414
## 6 2015 809. 1205. 0.672
## 7 2016 781. 1435. 0.544
## 8 2017 887. 1514. 0.586
## 9 2018 807. 1314. 0.614
## 10 2019 811. 1512. 0.537
## 11 2020 832. 1445. 0.576
## 12 2021 819. 1403. 0.584
## 13 2022 862. 1151. 0.749
# Plot evapotranpitation to precipitation ratio
plot(evapotrans_prec$Year, evapotrans_prec$Ratio, type = "b", col = "orange",
xlab = "Year", ylab = "Ratio", main = "Evapotranspiration to Precipitation Ratio")
The evapotranspiration to precipitation ratio demonstrates annual variations, ranging from approximately 0.41 to 0.75. Years like 2011 and 2022 show higher ratios, indicating a larger fraction of precipitation is lost to evapotranspiration. Conversely, years like 2014 and 2010 have lower ratios, suggesting more precipitation is available for runoff and infiltration. These fluctuations highlight the influence of climatic variability on water balance dynamics within the catchment.
This section focuses
on low-flow and trend analysis of the hydrological data. We will conduct
a low-flow analysis using the lfstat package, which
includes baseflow separation and plotting a flow duration curve.
Additionally, we will calculate key low-flow indexes such as BFI
(Baseflow Index), MAM (Mean Annual Minimum flow), and quantile flows
(Q95, Q90, Q70). The section will also explore drought periods based on
the Q95 threshold and analyze seasonal trends, including the
Mann-Kendall trend test for variables like discharge, precipitation, air
temperature, and evapotranspiration.
# Convert the Date column to POSIXct format for time series handling
main_data$Date <- as.POSIXct(strptime(main_data$Date, format = "%Y-%m-%d"))
# Define a zoo object for Discharge data with Date as the index
Qzoo <- zoo(main_data$Discharge, main_data$Date)
# Convert zoo object to lfobj using createlfobj
discharge_timeseries <- createlfobj(ts(Qzoo), startdate = "01/01/2010", hyearstart = 1)
# Define the units for the lfobj data
setlfunit("m^3/s")
# Conduct baseflow separation
discharge_timeseries$baseflow <- baseflow(discharge_timeseries$flow)
# Display the first few rows of the discharge_timeseries object
head(discharge_timeseries)
## day month year flow hyear baseflow
## 1 1 1 2010 59.058 2010 NA
## 2 2 1 2010 53.925 2010 NA
## 3 3 1 2010 47.681 2010 NA
## 4 4 1 2010 41.040 2010 NA
## 5 5 1 2010 38.694 2010 NA
## 6 6 1 2010 36.966 2010 NA
plot(discharge_timeseries$flow, type = "l", ylab = "Flow")
lines(discharge_timeseries$baseflow, col = "green")
The low-flow analysis plot displays the discharge over time (black line) along with the baseflow (green line) derived from the baseflow separation process. A noticeable decreasing trend in total flow is evident over the years, which may indicate long-term changes in climatic conditions, reduced precipitation, or increased water abstraction. The baseflow component remains relatively stable, indicating the consistent contribution of groundwater to the streamflow, even as overall flow declines.
bfplot(discharge_timeseries)
## NULL
# Plot the Flow Duration Curve
fdc(discharge_timeseries)
## 100% 99% 98% 97% 96% 95% 94%
## FDC-quantiles 859.744 265.1749 198.5197 165.5311 139.5318 124.5457 111.4355
## 93% 92% 91% 90% 89% 88% 87%
## FDC-quantiles 101.1971 94.24816 87.01376 81.7343 75.5104 71.20372 67.21979
## 86% 85% 84% 83% 82% 81% 80%
## FDC-quantiles 63.0349 59.76245 57.06852 54.55028 52.1816 50.30966 47.915
## 79% 78% 77% 76% 75% 74% 73%
## FDC-quantiles 46.00891 44.48692 42.86936 41.7912 40.2515 39.03082 37.94931
## 72% 71% 70% 69% 68% 67% 66%
## FDC-quantiles 36.84004 36.01759 35.1008 34.39782 33.47836 32.65884 31.92196
## 65% 64% 63% 62% 61% 60% 59%
## FDC-quantiles 30.96445 30.20936 29.56432 28.97764 28.25367 27.558 26.99014
## 58% 57% 56% 55% 54% 53% 52%
## FDC-quantiles 26.24132 25.62853 25.02784 24.5076 24.01942 23.52938 22.99832
## 51% 50% 49% 48% 47% 46% 45%
## FDC-quantiles 22.40255 22.0315 21.56409 21.09428 20.65336 20.27288 19.79815
## 44% 43% 42% 41% 40% 39% 38%
## FDC-quantiles 19.38068 18.96778 18.5832 18.25586 17.8788 17.58499 17.24476
## 37% 36% 35% 34% 33% 32% 31%
## FDC-quantiles 16.95739 16.73312 16.48505 16.18994 15.93157 15.70604 15.4277
## 30% 29% 28% 27% 26% 25% 24%
## FDC-quantiles 15.1942 14.95626 14.73848 14.46052 14.20366 13.96625 13.71984
## 23% 22% 21% 20% 19% 18% 17%
## FDC-quantiles 13.48862 13.22834 12.97761 12.7024 12.50939 12.278 12.051
## 16% 15% 14% 13% 12% 11% 10%
## FDC-quantiles 11.78856 11.53025 11.31492 11.05477 10.793 10.61953 10.3591
## 9% 8% 7% 6% 5% 4% 3% 2%
## FDC-quantiles 10.027 9.69504 9.34229 8.95922 8.63805 8.29416 7.883 7.42382
## 1% 0%
## FDC-quantiles 6.87029 5.515
abline(h=Q95(discharge_timeseries), col = "green", lty = 2)
text(90, 10, "Q95")
The Flow Duration Curve (FDC) shows the distribution of flow magnitudes over time in the catchment. The Q95 threshold (green dashed line) represents the flow that is exceeded 95% of the time, indicating low-flow conditions. The steep slope of the FDC at higher flows indicates a high variability of peak flows, while the flatter curve at lower flows reflects a stable baseflow regime. This result demonstrates the dominance of low to moderate flows in the catchment, with occasional high-flow events contributing to the variability observed in the upper range of the curve.
# FDC for the year 2014
y2014 <- fdc(discharge_timeseries, year = 2014)
# FDC for the year 2022
y2022 <- fdc(discharge_timeseries, year = 2022)
# Plot Flow Duration Curves (FDC) for 2014 and 2022
plot(as.numeric(y2014), type = "l", col = "blue", ylab = "Flow", main = "Flow Frequency Comparison: 2014 and 2022")
lines(as.numeric(y2022), col = "red", type = "l", lty = 2)
legend("topright", legend = c("2014", "2022"), col = c("blue", "red"), lty = c(1, 2))
The Flow Duration Curve (FDC) comparison between 2014 (blue line) and 2022 (red dashed line) shows a clear decrease in flow magnitudes over time. The FDC for 2014 indicates higher flows across all exceedance frequencies, suggesting wetter conditions and greater runoff compared to 2022. The lower flows in 2022 reflect a drier catchment, potentially due to reduced precipitation, increased water abstraction, or changes in land use. This comparison shows the progressive drying of the catchment and the variability in hydrological behavior over time.
# Calculate Baseflow Index
bfi <- BFI(discharge_timeseries)
bfi
## [1] 0.4925945
# Calculate Mean Annual Minimum Flow
mam <- MAM(discharge_timeseries)
mam
## [1] 8.581835
# Calculate Mean Flow
meanflow <- mean(discharge_timeseries$flow, na.rm = TRUE)
meanflow
## [1] 38.48014
# Calculate Flow Quantiles
# Q95
q95 <- Q95(discharge_timeseries)
q95
## Q95
## 8.63805
# Q90
q90 <- Q90(discharge_timeseries)
q90
## Q90
## 10.3591
# Q70
q70 <- Q70(discharge_timeseries)
q70
## Q70
## 15.1942
# Calculate Seasonality Index at Q95
seasonal_index <- seasindex(discharge_timeseries, Q = 95)
seasonal_index
## $theta
## [1] 3.857201
##
## $D
## [1] 224.0708
##
## $r
## [1] 0.7469678
# Calculate Seasonality Ratio
seasonal_ratio <- seasratio(discharge_timeseries, breakdays = "01/07", Q = 95)
seasonal_ratio
## Q95(1/1-01/07)/Q95(01/07-1/1)
## 1.3899
# Summarize
summary(discharge_timeseries)
## Startdate: 2010-01-01 (calendar year)
## Enddate: 2022-12-31 (calendar year)
##
## Time series covers 13 years.
## The hydrological year is set to start on January 1st.
##
## Meanflow MAM7 Q95 BFI
## 38.4800 8.5820 8.6380 0.4926
# Identify Drought Periods
discharge_timeseries$drought <- discharge_timeseries$flow < q95
# Display few rows of identified drought periods
drought_periods = discharge_timeseries[discharge_timeseries$drought == TRUE, ]
head(drought_periods)
## day month year flow hyear baseflow drought
## 191 10 7 2010 8.472 2010 8.4490 TRUE
## 192 11 7 2010 8.172 2010 8.1720 TRUE
## 193 12 7 2010 7.909 2010 7.9090 TRUE
## 194 13 7 2010 8.087 2010 7.8908 TRUE
## 195 14 7 2010 8.161 2010 7.8726 TRUE
## 196 15 7 2010 7.961 2010 7.8544 TRUE
# Create a Date column
drought_periods$Date <- as.Date(with(drought_periods, paste(year, month, day, sep = "-")), "%Y-%m-%d")
plot(drought_periods$Date, drought_periods$flow,
col = "#1a24a4", xlab = "Date", ylab = "Flow", main = "Flow with Drought Periods (Below Q95 Threshold)")
The drought periods are identified based on the Q95 threshold (8.638 m³/s), representing flows below which drought conditions occur. The plot shows variability in drought occurrences across the years. Drought events are minimal in 2010 and 2011 but increase significantly in 2012 and 2013, with prolonged low-flow conditions. No droughts are observed between 2014 and 2016 or 2019 and 2021, indicating periods of higher flows. In contrast, 2017 and 2022 show moderate to high numbers of drought days, reflecting a return to frequent low-flow conditions.
# Identify drought periods using the Q95 threshold
drought_time = find_droughts(discharge_timeseries, threshold = Q95(discharge_timeseries), na.rm = TRUE)
## Warning in as.xts.lfobj(x): No unit found in attributes, assuming 'm³/s'.
## Use flowunit(x) <- "l/s" to define the flow unit. See help(flowunit).
# Summarize drought events, filtering out minor ones
selected = summary(drought_time, drop_minor = 0)
# Plot the drought time series
plot(drought_time)
The plot displays drought periods identified using the Q95 threshold, with shaded segments indicating periods of low flow. The results show shorter drought durations around the end of 2010, followed by multiple droughts between 2012 and 2014, some of which are more prolonged. A significant drought period is observed toward the end of 2017, with even longer durations in the year 2022.
# Plot the seasonal distribution of selected drought events
plot(season(selected$time))
The plot illustrates the seasonal distribution of drought events. The results indicate that droughts predominantly occur during the summer season, followed by autumn. Spring experiences fewer droughts, while winter has the least occurrences due to lower evapotranspiration and stable baseflow conditions during colder months.
# Initialize list to store the results
mann_kandell_test <- list()
# Mann-Kendall test for Discharge
mann_kandell_test$Discharge <- MannKendall(main_data$Discharge)
# Mann-Kendall test for Precipitation at each station
precip_stations <- c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")
for (station in precip_stations) {
mann_kandell_test[[station]] <- MannKendall(main_data[[station]])
}
# Mann-Kendall test for Air Temperature
mann_kandell_test$Airtemp <- MannKendall(main_data$Airtemp)
# Mann-Kendall test for Evapotranspiration
mann_kandell_test$Evapotrans <- MannKendall(main_data$Evapotrans)
# Test result
mann_kandell_test
## $Discharge
## tau = -0.0314, 2-sided pvalue =0.0011597
##
## $P1_Luce
## tau = -0.0207, 2-sided pvalue =0.051273
##
## $P2_Solcava
## tau = -0.00822, 2-sided pvalue =0.44307
##
## $P3_Lasko
## tau = 0.000218, 2-sided pvalue =0.98394
##
## $P4_GornjiGrad
## tau = -0.0113, 2-sided pvalue =0.28991
##
## $P5_Radegunda
## tau = -0.00613, 2-sided pvalue =0.56852
##
## $Airtemp
## tau = 0.0207, 2-sided pvalue =0.032784
##
## $Evapotrans
## tau = 0.00605, 2-sided pvalue =0.53637
The Mann-Kendall test results indicate no significant trends for most precipitation stations, with P1_Luce showing a marginally decreasing trend. For air temperature, a slight increasing trend is observed with a statistically significant p-value, suggesting a gradual warming over time. Evapotranspiration shows no significant trend, indicating stable evapotranspiration levels throughout the analyzed period.
mk.test(main_data$Discharge)
##
## Mann-Kendall trend test
##
## data: main_data$Discharge
## z = -3.2486, n = 4748, p-value = 0.00116
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -3.543340e+05 1.189671e+10 -3.144278e-02
The Mann-Kendall test for discharge indicates a statistically significant decreasing trend, with a p-value of 0.00116. This result reflects a decline in water flow over time, suggesting potential long-term changes in hydrological conditions or water availability in the catchment.
plot(main_data$Date, main_data$Discharge, type = "l", xlab = "Year", ylab = "Discharge[m3/s]", main = "Discharge plot")
abline(lm(main_data$Discharge ~ main_data$Date), lty = 2, col = "red")
In this section, we
will perform a flood frequency analysis based on daily discharge data to
assess extreme flood events. The analysis begins with defining an annual
maxima sample derived from daily discharge records. Next, we will import
the vQvk peak discharge data, covering the complete
available dataset, and compare it with the annual maxima sample to
validate consistency. Using the vQvk data, we will conduct
a flood frequency analysis applying at least three different
distribution functions, and the results will be visualized, highlighting
measured discharge data under the Weibull distribution. Additionally,
confidence intervals will be calculated by the genci.simple
function from the lmomco package.
# Max discharge values
DischargeMax <- main_data %>%
group_by(Year) %>%
summarise(
maxDischarge = max(Discharge)
)
DischargeMax
## # A tibble: 13 × 2
## Year maxDischarge
## <dbl> <dbl>
## 1 2010 860.
## 2 2011 323.
## 3 2012 579.
## 4 2013 443.
## 5 2014 571.
## 6 2015 397.
## 7 2016 402.
## 8 2017 582.
## 9 2018 282.
## 10 2019 340.
## 11 2020 322.
## 12 2021 243.
## 13 2022 330.
# Plot yearly maximum discharge
plot(DischargeMax$Year, DischargeMax$maxDischarge, type = "l", xlab = "Year", ylab = "Max Discharge[m3/s]",
main = "Yearly Max Discharge")
abline(lm(DischargeMax$maxDischarge ~ DischargeMax$Year), lty = 2, col = "red")
The plot of yearly maximum discharge shows a clear declining trend in peak discharge values over the analyzed period. This is supported by the red dashed line, which represents the linear regression fit, further confirming the decreasing trend. The reduction in maximum discharge values could indicate a potential change in hydrological conditions, such as reduced precipitation intensity or changes in land use and catchment characteristics, affecting runoff generation during extreme events.
peak_data <- read_excel("peak_discharge_data.xlsx")
# Calculate annual peak for vQvk discharge
annual_maxima_vQvk <- peak_data %>%
group_by(Year) %>%
summarise(
vQvk = max(Qvk, na.rm = TRUE)
)
# Display few rows
head(annual_maxima_vQvk)
## # A tibble: 6 × 2
## Year vQvk
## <chr> <dbl>
## 1 1955 421
## 2 1956 569
## 3 1957 290
## 4 1958 524
## 5 1959 484
## 6 1960 502
# Filter vQvk data from 2010 onward
annual_vQvk_selected <- annual_maxima_vQvk %>%
filter(Year >= 2010)
# Combine the two peak discharge
combined_data <- data.frame(
Year = annual_vQvk_selected$Year,
main_DischargeMax = DischargeMax$maxDischarge,
vQvk_Discharge = annual_vQvk_selected$vQvk
)
# Determine y-axis limits
y_range <- range(c(combined_data$main_DischargeMax, combined_data$vQvk_Discharge), na.rm = TRUE)
# Plotting vQvk peak discharge
plot(combined_data$Year, combined_data$vQvk_Discharge, type = "l", col = "blue", lwd = 1,
xlab = "Year", ylab = "Peak Discharge", main = "Comparison of Peak Discharge Values", ylim = y_range)
# Add main_DischargeMax to the plot
lines(combined_data$Year, combined_data$main_DischargeMax, col = "red", lwd = 1)
# Adding a legend
legend("topright", legend = c("Main DischargeMax", "vQvk Discharge"),
col = c("red", "blue"), lwd = 1)
# Summary statistics for vQvk
summary(annual_maxima_vQvk)
## Year vQvk
## Length:66 Min. : 262.3
## Class :character 1st Qu.: 461.7
## Mode :character Median : 565.0
## Mean : 644.3
## 3rd Qu.: 766.8
## Max. :1490.0
# Plot annual maxima series for vQvk
plot(annual_maxima_vQvk$Year, annual_maxima_vQvk$vQvk, type = "b",
xlab = "Year", ylab = "vQvk [m3/s]", main = "Annual Peak Discharge for Savinja River")
# Calculate L-moments for the vQvk data
lmoments <- lmoms(annual_maxima_vQvk$vQvk)
# Display calculated L-moments
lmoments$lambdas
## [1] 644.270955 146.790664 41.173791 27.376282 5.118641
# Estimate parameters for PE3, GEV, GNO, Normal, and Weibull distributions
pe3par <- lmom2par(lmoments, type = "pe3")
gevpar <- lmom2par(lmoments, type = "gev")
gnopar <- lmom2par(lmoments, type = "gno")
norpar <- lmom2par(lmoments, type = "nor")
wei_par <- lmom2par(lmoments, type = "wei")
# Define return periods and calculate exceedance probabilities
return_periods <- c(2, 5, 10, 20, 30, 50, 100, 300,500)
# Calculate Fx values for each return period
exceedance_probs <- 1 - 1 / return_periods
# Calculate discharge values for each return period using PE3, GEV, GNO, Normal, and Weibull distributions
discharge_pe3 <- quape3(exceedance_probs, pe3par)
discharge_gev <- quagev(exceedance_probs, gevpar)
discharge_gno <- quagno(exceedance_probs, gnopar)
discharge_nor <- quanor(exceedance_probs, norpar)
discharge_wei <- quawei(exceedance_probs, wei_par)
# Calculate Weibull plotting positions for observed data and return periods for measured data
weibull_positions <- pp(annual_maxima_vQvk$vQvk, a = 0)
observed_return_periods <- 1 / (1 - weibull_positions)
# Plot the PE3, GEV, GNO, Normal, and Weibull distribution results
plot(return_periods, discharge_pe3, type = "l", log = "x", col = "blue",
xlab = "Return Period [years]", ylab = "Discharge [m3/s]",
main = "Flood Frequency Analysis for Savinja River", ylim = c(400, max(discharge_pe3, discharge_gev, discharge_gno, discharge_wei, na.rm = TRUE)))
# Add GEV, GNO, Normal, and Weibull distribution lines
lines(return_periods, discharge_gev, col = "red")
lines(return_periods, discharge_gno, col = "green")
lines(return_periods, discharge_nor, col = "orange")
lines(return_periods, discharge_wei, col = "#ff31a1")
# Add measured discharge data points
points(observed_return_periods, sort(annual_maxima_vQvk$vQvk), col = "black", pch = 16)
# Add legend
legend("topleft", legend = c("PE3 Distribution", "GEV Distribution", "GNO Distribution",
"Normal Distribution", "Weibull Distribution", "Observed Data"),
col = c("blue", "red", "green", "orange", "#ff31a1", "black"), lty = 1, pch = c(NA, NA, NA, NA, NA))
The flood frequency analysis illustrates discharge values for various return periods using multiple probability distributions, including PE3, GEV, GNO, Normal, and Weibull. The Weibull and PE3 distributions show a closer fit to the observed data, particularly for low and moderate return periods, making them more reliable for this dataset. On the other hand, GEV and GNO distributions tend to overestimate discharge values, especially for higher return periods, while the Normal distribution consistently underestimates discharge across most return periods. This comparison shows the importance of selecting an appropriate distribution for accurate flood frequency analysis in hydrological studies.
# Set seed for reproducibility
set.seed(123)
# Calculate confidence intervals for the PE3 distribution
ci_pe3 <- genci.simple(para = pe3par, n = length(annual_maxima_vQvk$vQvk),
f = exceedance_probs, level = 0.90, nsim = 1000, edist = "gno")
## 9-8-7-6-5-4-3-2-1-0-
# Combine the confidence intervals and design discharges
confidence_intervals <- data.frame(
Return_Period = return_periods,
Lower_CI = ci_pe3$lwr,
Upper_CI = ci_pe3$upr,
PE3_Discharge = discharge_pe3
)
# Plot PE3 distribution with confidence intervals
plot(confidence_intervals$Return_Period, confidence_intervals$PE3_Discharge, type = "l", log = "x",
col = "green", xlab = "Return Period [years]", ylab = "Discharge [m3/s]",
main = "Flood Frequency Analysis",
ylim = c(400, max(ci_pe3$upr, na.rm = TRUE)))
# Add confidence interval bounds for PE3 distribution
lines(confidence_intervals$Return_Period, confidence_intervals$Lower_CI, col = "purple", lty = 2)
lines(confidence_intervals$Return_Period, confidence_intervals$Upper_CI, col = "purple", lty = 2)
# Add observed discharge data points for comparison
points(observed_return_periods, sort(annual_maxima_vQvk$vQvk), col = "black", pch = 16)
# Add legend
legend("bottomright", legend = c("PE3 Distribution", "Lower CI", "Upper CI", "Observed Data"),
col = c("green", "purple", "purple", "black"), lty = c(1, 2, 2, NA), pch = c(NA, NA, NA, 16))
Based on the previous analysis, the Weibull and PE3 distributions were observed to fit the dataset well, with PE3 being slightly more consistent with the observed data. As a result, the PE3 distribution was selected for calculating the confidence intervals. The plot shows the PE3 discharge curve (green line) along with the 90% confidence intervals (purple dashed lines). The observed data points (black dots) are mostly contained within the confidence intervals, demonstrating the reliability of the PE3 model. Additionally, the confidence intervals widen for larger return periods, reflecting increased uncertainty in predicting extreme events such as rare floods.
In this section, we
will analyze precipitation data to explore trends and variability across
multiple thresholds and indices. We begin by calculating the total
number of days with measurable precipitation, as well as days with
precipitation exceeding 10 mm and 50 mm, followed by determining the
probability of occurrence for each threshold. Next, we will compute
Standardized Precipitation Index (SPI) values over 1, 3, 6, and 12-month
periods for each of the five precipitation stations to assess drought
conditions. The identified drought periods will then be compared to
those derived from discharge data using the lfstat package,
providing insight into the relationship between precipitation and
drought events. Additionally, we will calculate the Effective Drought
(ED) index and examine its relationship with the SPI results. Finally,
we will incorporate the North Atlantic Oscillation (NAO) index by
downloading NAO data and calculating its correlation with precipitation
records to investigate potential climate influences on regional
precipitation patterns.
# Calculate days with precipitation > 10 mm
prep_above_10 <- main_data %>%
summarise(
P1_Luce = sum(P1_Luce > 10),
P2_Solcava = sum(P2_Solcava > 10),
P3_Lasko = sum(P3_Lasko > 10),
P4_GornjiGrad = sum(P4_GornjiGrad > 10),
P5_Radegunda = sum(P5_Radegunda > 10),
)
# Calculate frequency for days with precipitation > 10 mm
frequency_above_10 <- prep_above_10 * 100 / nrow(main_data)
# Combine both the counts and frequency rows
results_above_10 <- bind_rows(prep_above_10, frequency_above_10)
rownames(results_above_10) <- c("Days", "Frequency")
# Display the results
results_above_10
## P1_Luce P2_Solcava P3_Lasko P4_GornjiGrad P5_Radegunda
## Days 624.00000 627.00000 538.00000 640.00000 589.00000
## Frequency 13.14238 13.20556 11.33109 13.47936 12.40522
# Calculate days with precipitation > 50 mm
prep_above_50 <- main_data %>%
summarise(
P1_Luce = sum(P1_Luce > 50),
P2_Solcava = sum(P2_Solcava > 50),
P3_Lasko = sum(P3_Lasko > 50),
P4_GornjiGrad = sum(P4_GornjiGrad > 50),
P5_Radegunda = sum(P5_Radegunda > 50)
)
# Calculate frequency for days with precipitation > 50 mm
frequency_above_50 <- prep_above_50 * 100 / nrow(main_data)
# Combine both the counts and frequency rows
results_above_50 <- bind_rows(prep_above_50, frequency_above_50)
rownames(results_above_50) <- c("Days", "Frequency")
# Display the results
results_above_50
## P1_Luce P2_Solcava P3_Lasko P4_GornjiGrad P5_Radegunda
## Days 59.000000 62.000000 18.000000 50.000000 31.0000000
## Frequency 1.242628 1.305813 0.379107 1.053075 0.6529065
# Create a zoo object for all 5 stations
precip_zoo <- zoo(main_data[, c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda")], order.by = main_data$Date)
head(precip_zoo)
## P1_Luce P2_Solcava P3_Lasko P4_GornjiGrad P5_Radegunda
## 2010-01-01 7.0 4.8 4.3 4.6 2.2
## 2010-01-02 3.0 0.6 2.1 3.2 2.0
## 2010-01-03 0.1 2.7 1.5 0.0 0.1
## 2010-01-04 0.0 0.0 0.0 0.0 0.0
## 2010-01-05 3.3 2.0 7.0 3.1 4.0
## 2010-01-06 6.2 4.6 7.0 9.2 8.2
# Convert time index to year-month
yearmon_index <- as.yearmon(time(precip_zoo))
# Aggregate precipitation data by month
monthly_precip <- aggregate(precip_zoo, yearmon_index, sum)
# Plot the aggregated monthly precipitation data for P1_Luce
plot(monthly_precip$P1_Luce, xlab="Year", ylab="Monthly precipitation [mm]", main="Monthly Precipitation P1_Luce")
# Get summary of the monthly data
summary(monthly_precip)
## Index P1_Luce P2_Solcava P3_Lasko
## Min. :2010 Min. : 0.1 Min. : 0.3 Min. : 0.00
## 1st Qu.:2013 1st Qu.: 66.5 1st Qu.: 61.0 1st Qu.: 55.90
## Median :2016 Median :117.3 Median :112.4 Median : 87.50
## Mean :2016 Mean :128.1 Mean :128.0 Mean : 98.45
## 3rd Qu.:2020 3rd Qu.:174.7 3rd Qu.:185.1 3rd Qu.:135.80
## Max. :2023 Max. :382.7 Max. :398.0 Max. :313.20
## P4_GornjiGrad P5_Radegunda
## Min. : 0.0 Min. : 0.0
## 1st Qu.: 67.2 1st Qu.: 63.6
## Median :113.4 Median :103.0
## Mean :124.6 Mean :114.5
## 3rd Qu.:166.6 3rd Qu.:166.4
## Max. :360.6 Max. :326.0
# Calculate SPI at 3-month scale for all 5 stations
# P1_Luce SPI (3-month scale)
spi_3_P1_Luce <- spi(data = as.numeric(monthly_precip$P1_Luce), scale = 3, distribution = "Gamma")
## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."
# P2_Solcava SPI (3-month scale)
spi_3_P2_Solcava <- spi(data = as.numeric(monthly_precip$P2_Solcava), scale = 3, distribution = "Gamma")
## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."
# P3_Lasko SPI (3-month scale)
spi_3_P3_Lasko <- spi(data = as.numeric(monthly_precip$P3_Lasko), scale = 3, distribution = "Gamma")
## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."
# P4_GornjiGrad SPI (3-month scale)
spi_3_P4_GornjiGrad <- spi(data = as.numeric(monthly_precip$P4_GornjiGrad), scale = 3, distribution = "Gamma")
## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."
# P5_Radegunda SPI (3-month scale)
spi_3_P5_Radegunda <- spi(data = as.numeric(monthly_precip$P5_Radegunda), scale = 3, distribution = "Gamma")
## [1] "Calculating the Standardized Precipitation Evapotranspiration Index (SPEI) at a time scale of 3. Using kernel type 'rectangular', with 0 shift. Fitting the data to a Gamma distribution. Using the ub-pwm parameter fitting method. Checking for missing values (`NA`): all the data must be complete. Using the whole time series as reference period. Input type is vector. No time information provided, assuming a monthly time series."
# Plot SPI-3 for P1_Luce
plot(spi_3_P1_Luce, main = "SPI-3 for P1_Luce")
## Warning: Removed 155 rows containing missing values (`geom_point()`).
# Plot SPI-3 for P2_Solcava
plot(spi_3_P2_Solcava, main = "SPI-3 for P2_Solcava")
## Warning: Removed 155 rows containing missing values (`geom_point()`).
# Plot SPI-3 for P3_Lasko
plot(spi_3_P3_Lasko, main = "SPI-3 for P3_Lasko")
## Warning: Removed 155 rows containing missing values (`geom_point()`).
# Plot SPI-3 for P4_GornjiGrad
plot(spi_3_P4_GornjiGrad, main = "SPI-3 for P4_GornjiGrad")
## Warning: Removed 155 rows containing missing values (`geom_point()`).
# Plot SPI-3 for P5_Radegunda
plot(spi_3_P5_Radegunda, main = "SPI-3 for P5_Radegunda")
## Warning: Removed 155 rows containing missing values (`geom_point()`).
# View SPI values
spi_3_P1_Luce
## [1] NA NA -0.7204634761 -0.1035448200 -0.5539357550
## [6] -0.6676236843 -1.2969851317 -0.0655589479 0.3259998529 1.7653594567
## [11] 1.1743821379 1.0562537491 0.7864004975 0.5315218041 0.0836068510
## [16] -0.3778538826 -0.1503402766 0.0434250727 0.0874669662 0.5337545988
## [21] 0.4556579062 -0.4739078884 0.1009485414 -1.0554242865 -0.5487552819
## [26] -1.9509528509 -0.7768173785 -1.8320052570 -1.2106029338 0.0593617193
## [31] 0.7865403272 0.0111125757 0.4313171070 1.0155799377 1.8984148098
## [36] 2.0180754590 1.6783405425 0.3900624293 0.2186045794 0.9304389929
## [41] 1.1681300011 1.5943340047 -0.2344968859 -0.7152876314 -2.0219231681
## [46] -1.6088188273 -0.8308576836 0.2313359447 0.3006416530 1.4541264415
## [51] 1.8924172585 1.7712864233 1.3287712361 -0.9517816819 -0.1692218462
## [56] 0.4717715631 1.7397842350 1.5117932583 0.7086622877 0.8132757906
## [61] 0.5325371731 0.5917726450 -0.6862657627 -0.1845212838 -1.0187620935
## [66] -0.5429304619 -0.2155548230 1.1056130145 0.4380425893 0.0857337977
## [71] -0.0005738679 -0.5370730805 -2.0228457110 -1.2656144437 0.7686494723
## [76] 1.2674089475 1.0134699039 -0.6940074221 -0.6938430762 -0.0260690707
## [81] 0.6011259267 -0.7656586966 -0.8326843161 -0.1094490102 -0.0634693010
## [86] -0.2205900906 -0.8041570121 -0.3723307295 1.1224227183 0.0385131016
## [91] 1.0044809018 -1.7667037986 -0.6991236156 0.3616881567 0.1583698010
## [96] 0.7313631578 0.8082305226 0.9764092742 0.9127374322 0.5074158067
## [101] 0.7981022991 1.0834361643 0.7710131441 -0.0281412397 0.0243632898
## [106] -0.6675149485 0.2305222973 -0.5486176001 -1.1294044301 -1.4507829219
## [111] -0.6551449107 -0.2509881927 0.0996966452 0.3182928461 0.1555393112
## [116] 0.1787635734 -0.2419443280 -0.2190713667 -1.0221250319 0.1876205391
## [121] 0.7997899101 0.6987957449 -0.2889669808 -0.3224278584 -0.6310426306
## [126] -0.4581669440 -0.0433886664 0.7751208418 1.0284666199 0.0982773229
## [131] 0.5755299204 -0.7070806847 0.3580075641 0.4135041528 1.2872524581
## [136] 0.4741089368 -0.6088379641 1.6954463058 1.8524733713 1.6152503773
## [141] -0.4130392476 -0.6971142888 -1.5138998942 -1.2742193808 -0.6216554028
## [146] -0.1473223838 -0.7668758330 -1.2083643412 -1.0946375692 -1.3541367833
## [151] -1.9365516846 -2.0665831884 -1.6110726896 -0.2567218298 -0.4366623883
## [156] -0.6084039781 -0.7632183398
spi_3_P2_Solcava
## [1] NA NA -0.709959804 -0.174533764 -0.497014226
## [6] -0.430933635 0.011703359 -0.002330165 0.482201164 0.897560815
## [11] 1.009847864 0.855183314 0.678782311 0.335845732 -0.121475490
## [16] -0.440607366 -0.323648029 0.215482312 0.756181542 1.017861769
## [21] 0.507424063 -0.170113402 0.073484103 -0.986796098 -1.101574997
## [26] -1.685578121 -0.840649208 -1.767250300 -1.317421508 -0.280943819
## [31] 0.126513193 -0.238080607 0.149988713 0.986911951 1.645075609
## [36] 1.783900628 1.110181653 0.172371864 0.028718380 0.805148911
## [41] 0.928405189 1.732601397 0.030641832 -0.198973539 -1.700147450
## [46] -1.034814771 -0.666018796 0.580198655 0.757928699 1.443396110
## [51] 1.881046606 1.968251405 1.358931228 -0.333713907 1.044196045
## [56] 1.090161763 1.702052692 1.263730411 0.692481686 0.985183101
## [61] 0.388022824 0.574388943 -0.416976499 0.132888314 -0.609506906
## [66] 0.035820221 0.344004865 1.160867857 0.251307581 0.145834550
## [71] 0.264800164 -0.276429035 -1.785874198 -1.499691736 0.385879799
## [76] 1.057238500 1.057162577 -0.528289628 -0.771281198 -0.760759640
## [81] -0.095289536 -1.464799964 -1.152047817 -0.817679789 -0.688842188
## [86] -0.467649344 -1.023298834 -0.584115626 0.895222854 -0.355491530
## [91] 0.216427568 -2.139273378 -0.727013158 0.413397043 0.173568293
## [96] 0.575557501 0.833935448 0.918890410 0.980835666 0.327726682
## [101] 0.496487980 0.467792171 -0.200755622 -0.684921188 -0.292476496
## [106] -1.181152146 -0.072779909 -0.955516915 -1.005298359 -1.203718306
## [111] -0.194040103 0.234508478 1.031230752 1.618783720 1.495363490
## [116] 0.312359304 -0.842060763 -0.577345665 -1.165186252 0.709291848
## [121] 1.214961313 0.960651077 -0.273540815 -0.340352115 -0.684128580
## [126] -1.096569331 -0.957714796 0.863493787 1.727339804 1.220022235
## [131] 1.170118347 -0.278758781 0.872198149 0.633472779 1.471187346
## [136] 0.330420936 -0.835875828 0.753879399 0.686983065 0.778310548
## [141] -0.555754956 -0.818326861 -1.264092922 -1.237272905 -0.510860093
## [146] -0.052460659 -0.604088821 -1.230098127 -1.273425239 -1.668011192
## [151] -2.823708556 -1.119114300 -0.462859516 0.495853942 -0.476376599
## [156] -0.743817312 -0.644799247
spi_3_P3_Lasko
## [1] NA NA -0.318730311 0.488293739 0.065516561
## [6] 0.140802544 -0.266518721 -0.148910456 -0.106848116 1.409528093
## [11] 1.510970123 1.456046294 0.766348515 0.638035797 -0.852248834
## [16] -0.829956754 -0.991936943 -0.828203542 -0.293253544 0.529311600
## [21] -0.031549226 -1.359556061 -1.257381728 -1.281806457 -0.595221145
## [26] -1.436736644 -0.521976870 -1.459844374 -0.725496543 -0.260814555
## [31] 0.330452767 -0.598891737 -0.585915291 0.031033203 1.496063092
## [36] 1.255116717 1.358655303 0.268482679 1.230483594 1.539408872
## [41] 1.527435927 1.431104275 -0.575285249 -1.580625020 -2.284673982
## [46] -1.047000350 -0.570143401 0.842604951 0.500459791 1.491209131
## [51] 1.285053884 0.992755907 1.170389096 -0.120553861 1.158970685
## [56] 0.573431040 1.373631969 1.044575667 0.558551282 -0.109421466
## [61] -0.557813884 0.112793078 0.215496499 0.016989098 -0.703102756
## [66] -0.283211706 1.088779137 0.955346370 0.355488130 -1.111079187
## [71] -0.009017051 -0.249145716 -0.708948391 -1.433410298 0.831680932
## [76] 1.209465025 1.306399922 0.650537236 0.791066622 0.086881572
## [81] 0.207962502 -0.968355731 -0.613117465 -0.352262698 -0.362808855
## [86] -0.346840636 -1.010045864 -0.346657549 0.078569339 -1.251390684
## [91] -1.265874200 -2.300301324 -1.362828095 0.386930099 0.680867061
## [96] 1.221848345 1.168069474 1.393903874 1.343450141 0.952756451
## [101] 0.967233236 0.507392742 0.409086977 0.210255956 0.305460454
## [106] -0.833081099 -1.337074546 -1.693878254 -2.922012442 -1.415180762
## [111] -1.280614179 -0.199270202 -0.129750387 1.176577925 1.127680993
## [116] 0.890921681 0.717808484 0.480683148 -0.002944508 -0.198057651
## [121] 0.322843367 0.219324308 -0.808543915 -1.022561747 -1.261416217
## [126] -1.715441140 -2.355398715 0.795597006 1.561060027 1.146513734
## [131] 0.516280430 -0.393488565 0.792705714 0.305967880 0.573673561
## [136] -0.204531595 -0.490676795 1.174080579 0.332098985 1.253281814
## [141] 0.320224062 0.306539246 -1.056993646 -0.591662718 0.138702313
## [146] 0.321065691 -0.386549797 -0.686029187 -0.457444061 -0.446384696
## [151] -0.447952062 -0.659201058 -0.432661621 0.680414206 0.278063602
## [156] 0.286271663 0.001287518
spi_3_P4_GornjiGrad
## [1] NA NA -0.717472654 0.076834022 -0.509205904
## [6] -0.641802564 -1.137559357 -0.277074732 0.052024124 1.730623873
## [11] 1.401837588 1.465811714 0.918222709 0.629490871 -0.092858565
## [16] -0.487177260 -0.772731759 -0.330128553 -0.910541763 0.312566750
## [21] -0.085461014 -1.160102389 -0.434884080 -1.443658177 -0.660788358
## [26] -1.731770298 -0.665023775 -1.612549391 -1.052019598 -0.537461031
## [31] 0.100328925 -0.405067227 0.319167925 0.337925852 1.701775205
## [36] 1.819079628 1.886701377 0.430722486 0.307797916 0.945729254
## [41] 1.086669978 1.590199261 -0.419707846 -0.809905387 -1.668256578
## [46] -1.188379664 -0.709739296 0.399384613 0.086265904 1.419926217
## [51] 1.756156175 1.661769776 1.077692434 -1.173949001 0.350411446
## [56] 0.775397181 2.031604855 1.492024916 0.441760628 0.420299822
## [61] 0.170550693 0.561980492 -0.540429925 -0.384328244 -1.089042964
## [66] -0.273592115 -0.120395372 0.846604922 -0.359493840 -0.516441045
## [71] 0.017549446 -0.364089955 -1.546727260 -1.457216820 0.703219369
## [76] 1.213283386 1.239926852 0.203678309 0.448897549 0.694729658
## [81] 0.307810096 -1.101537535 -1.031462389 -0.293718858 -0.389226690
## [86] -0.507335114 -0.767202848 -0.256958543 1.037415325 -0.039477339
## [91] 1.597655885 -1.868749127 -0.459847253 0.456893742 0.189610133
## [96] 0.646753503 0.785376096 1.010224670 1.154860798 0.746574675
## [101] 1.177492736 1.453076062 0.302829083 0.286010087 0.199736572
## [106] -0.232931152 -0.111420133 -0.614940656 -1.170389230 -1.260628318
## [111] -0.738937044 -0.252535653 0.084382412 0.373797320 1.076504481
## [116] 0.854066143 0.184266943 -0.046778143 -1.033464104 0.089868710
## [121] 0.679846076 0.655953201 -0.257905162 -0.247386060 -0.536325046
## [126] -0.541770764 -0.804365514 0.558123086 1.386605466 0.716869443
## [131] 1.042993738 -0.481366704 0.541156220 0.526610814 1.285417945
## [136] 0.427331563 -0.376722118 1.366523963 1.156068197 1.170155741
## [141] -0.154997473 -0.549505026 -1.247701950 -1.327269842 -0.340562236
## [146] -0.207353722 -1.014767545 -1.507194995 -1.150485306 -1.286264408
## [151] -1.540854110 -2.148708436 -1.699026634 0.233382711 -0.004262021
## [156] -0.131077374 -0.758694226
spi_3_P5_Radegunda
## [1] NA NA -0.46053859 0.19529862 -0.56052002 -1.38783093
## [7] -2.18017336 -1.39502840 -0.23133538 1.43361604 1.41272441 1.25066023
## [13] 0.50004722 0.32821903 -0.34139795 -0.57117429 -0.66483260 -0.35792566
## [19] 0.01956341 0.37962939 0.16036838 -1.18372706 -0.51010024 -1.35839118
## [25] -0.71174862 -1.66791004 -0.77175665 -1.47716360 -1.12744349 -0.53879633
## [31] -0.34788656 -0.40451088 0.36065990 0.62432914 1.77088452 1.62800829
## [37] 1.79809306 0.47194107 0.65040389 1.38141828 1.64426005 1.65956725
## [43] -0.80097369 -1.55638375 -2.36227460 -1.12382972 -0.37509596 0.94221759
## [49] 0.76132537 1.60930134 1.77986025 1.42473820 1.02569029 -0.64434628
## [55] 0.46654481 0.57965064 1.77813259 1.47034757 0.75836076 0.59761154
## [61] 0.39920572 0.57927831 -0.52586587 -0.18738369 -0.61794857 0.12206524
## [67] 0.51922133 1.02863770 -0.03497377 -0.79368756 -0.30577124 -0.44947984
## [73] -1.36593051 -1.38925051 0.71747277 1.10143005 1.11323843 -0.29095495
## [79] 0.57756668 0.69182793 0.54327157 -1.37012381 -1.15834451 -0.65424824
## [85] -0.41275782 -0.69167089 -0.87962541 -0.58650300 0.62730506 -0.33372181
## [91] 1.81379217 -0.58206517 0.02491121 0.33043076 0.26588199 0.83927758
## [97] 0.46424837 0.88906381 1.04890927 0.86248803 1.15152911 1.27242044
## [103] 0.52072197 -0.37928760 0.11726215 -0.01987118 0.37430108 -0.83615355
## [109] -1.67353911 -1.22577520 -0.98860032 -0.23849151 -0.09396873 0.68324493
## [115] 0.34958355 0.83340617 0.12858089 0.26863449 -0.93984337 0.11954809
## [121] 0.81382205 0.70853189 -0.08534530 -0.35068384 -0.73831988 -0.83540773
## [127] -0.96464894 1.24462259 1.41200438 0.76549531 0.44217202 -0.45960842
## [133] 0.59028135 0.40284506 1.13664198 0.19383232 -0.65846001 1.44140306
## [139] 0.85236465 0.87667259 -0.92922779 -0.52770821 -1.33014400 -0.93796357
## [145] -0.19290923 0.09338470 -0.86253475 -1.33605039 -0.84822271 -0.62448379
## [151] -0.78182356 -1.22239232 -0.96133671 0.29628077 -0.20292438 -0.46606478
## [157] -0.82313394
# Convert data to daily for P1_Luce
daily_date <- as.Date(time(precip_zoo))
# Daily precipitation totals for P1_Luce
daily_precip <- aggregate(precip_zoo$P1_Luce, daily_date, sum)
# Check the basic statistics of the daily data
summary(daily_precip)
## Index daily_precip
## Min. :2009-12-31 Min. : 0.000
## 1st Qu.:2013-03-31 1st Qu.: 0.000
## Median :2016-06-30 Median : 0.000
## Mean :2016-06-30 Mean : 4.236
## 3rd Qu.:2019-09-30 3rd Qu.: 2.200
## Max. :2022-12-30 Max. :116.100
# Plot daily precipitation to examine data structure
plot(daily_precip, xlab="Year", ylab="Daily precipitation [mm]", main="Daily Precipitation for P1_Luce")
# Parameters
observation_days <- 365
cumulative_precip <- rep(NA, length(daily_precip))
# Convert daily precipitation data to a numeric vector
precip_values <- as.numeric(daily_precip)
# Calculate cumulative precipitation (EP)
for (day in (observation_days + 1):length(precip_values)) {
cumulative_sum <- 0
for (past_day in 1:observation_days) {
cumulative_sum <- cumulative_sum + sum(precip_values[(day - past_day + 1):day]) / past_day
}
cumulative_precip[day] <- cumulative_sum
}
# Mean of cumulative precipitation
mean_cumulative_precip <- mean(cumulative_precip, na.rm=TRUE)
# Deviation from mean
deviation_precip <- cumulative_precip - mean_cumulative_precip
# Standardized EDI
EDI <- deviation_precip / sd(deviation_precip, na.rm=TRUE)
# Plot EDI with drought and wetness severity lines
plot(time(daily_precip), EDI, type="l", col="#000000", lty=3,
xlab="Year", ylab="EDI", main="Effective Drought Index (EDI) for P1_Luce", xaxt="n")
# Define yearly labels for the x-axis
years <- seq(2010, 2022, by = 1)
x_positions <- seq(from = as.numeric(time(daily_precip)[1]),
to = as.numeric(time(daily_precip)[length(time(daily_precip))]),
length.out = length(years))
# Add x-axis with yearly labels
axis(1, at = x_positions, labels = years)
# Add threshold lines
abline(h = c(-2, -1.5, -1, 1, 1.5, 2),
col = c("red", "orange", "yellow", "green", "lightblue", "blue"),
lty = 2)
# Add labels for each threshold
text(x = rep(mean(x_positions), 6), y = c(-2, -1.5, -1, 1, 1.5, 2),
labels = c("Extremely Dry", "Severely Dry", "Moderately Dry",
"Moderately Wet", "Severely Wet", "Extremely Wet"),
col = c("red", "orange", "yellow", "green", "lightblue", "blue"),
pos = 4)
The plot shows the Effective Drought Index (EDI) values for P1_Luce, highlighting variations in drought and wetness conditions from 2010 to 2022. The years 2012 and mid-2021 display extremely dry conditions (EDI below -2), while mid-2012, late 2013, and mid-2014 indicate extremely wet conditions (EDI above 2). Seasonal fluctuations are evident, reflecting the influence of annual weather patterns.
# Plot SPI for comparison
par(mar = c(3, 3, 3, 3) + 0.1)
plot(time(monthly_precip), spi_3_P1_Luce$fitted, type="l", col="#000000", lty=3,
xlab="Year", ylab="SPI-3", main="Standardized Precipitation Index (SPI-3) for P1_Luce")
# Add threshold lines
abline(h = c(-2, -1.5, -1, 1, 1.5, 2),
col = c("red", "orange", "yellow", "green", "lightblue", "blue"),
lty = 2)
# Add labels for each threshold
text(x = rep(2011, 6), y = c(-2, -1.5, -1, 1, 1.5, 2),
labels = c("Extremely Dry", "Severely Dry", "Moderately Dry", "Moderately Wet", "Severely Wet", "Extremely Wet"),
col = c("red", "orange", "yellow", "green", "lightblue", "blue"))
drought_zoo = zoo(selected$duration, selected$time)
par(new = TRUE)
plot(drought_zoo, axes = FALSE, xlab = "", ylab = "", type = "h", col = "purple", lwd = 2)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a
## graphical parameter
axis(side = 4, at = pretty(range(drought_zoo)))
mtext("Duration of hydrological drought", side = 4, line = 3)
The plot compares the SPI-3 values for P1_Luce with the previously identified drought periods derived from the discharge data. Years such as 2010, 2011, 2012, 2017, and 2022 show consistency between meteorological droughts (negative SPI-3) and hydrological droughts (purple bars). However, some discrepancies are observed, such as in mid-2013, where precipitation data indicates wet conditions (positive SPI-3), while discharge data identifies drought conditions. This differences in meteorological and hydrological droughts are captured and emphasizes the complementary nature of these indices in understanding drought dynamics.
# Download NAO index data
nao <- download_nao()
# Check the structure of the data
head(nao)
## # A tibble: 6 × 3
## Year Month NAO
## <int> <ord> <dbl>
## 1 1950 Jan 0.92
## 2 1950 Feb 0.4
## 3 1950 Mar -0.36
## 4 1950 Apr 0.73
## 5 1950 May -0.59
## 6 1950 Jun -0.06
# Basic statistics
summary(nao$NAO)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.180000 -0.750000 0.040000 -0.004378 0.720000 3.040000
# Plot NAO index data
plot(nao$NAO,ylab="NAO",type="l")
# Autocorrelation in NAO data
acf(nao$NAO[1:200])
# Aggregate Precipitation Data by Month
monthly_precip <- main_data %>%
group_by(Year, Month) %>%
summarize(
P1_Luce = mean(P1_Luce, na.rm = TRUE),
P2_Solcava = mean(P2_Solcava, na.rm = TRUE),
P3_Lasko = mean(P3_Lasko, na.rm = TRUE),
P4_GornjiGrad = mean(P4_GornjiGrad, na.rm = TRUE),
P5_Radegunda = mean(P5_Radegunda, na.rm = TRUE)
) %>%
mutate(YearMonth = as.yearmon(paste(Year, Month), "%Y %m"))
## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.
# Prepare NAO Data
nao <- nao %>%
mutate(YearMonth = as.yearmon(paste(Year, Month), "%Y %b"))
# Merge Precipitation and NAO Data
merged_data <- merge(monthly_precip, nao, by = "YearMonth", all.x = TRUE)
# Compute Correlations
correlations <- merged_data %>%
dplyr::select(P1_Luce, P2_Solcava, P3_Lasko, P4_GornjiGrad, P5_Radegunda, NAO) %>%
dplyr::summarize(across(-NAO, ~ cor(., merged_data$NAO))) %>%
tidyr::pivot_longer(cols = everything(), names_to = "Station", values_to = "Correlation")
# Display correlations
correlations
## # A tibble: 5 × 2
## Station Correlation
## <chr> <dbl>
## 1 P1_Luce -0.231
## 2 P2_Solcava -0.228
## 3 P3_Lasko -0.286
## 4 P4_GornjiGrad -0.217
## 5 P5_Radegunda -0.257
# Visualize correlations
ggplot(correlations, aes(x = Station, y = Correlation, fill = Station)) +
geom_bar(stat = "identity") +
labs(
title = "Correlation between NAO and Precipitation",
x = "Precipitation Station",
y = "Correlation Coefficient"
)
The bar chart show the correlation coefficients between the NAO (North Atlantic Oscillation) index and precipitation data across five meteorological stations. All stations exhibit a negative correlation, with values ranging from -0.217 at P4_GornjiGrad to -0.286 at P3_Lasko. This suggests that an increase in the NAO index is generally associated with decreased precipitation at these stations. The magnitude of these correlations highlights the varying sensitivity of different stations to atmospheric conditions influenced by the NAO.
This section focuses on analyzing the spatial characteristics of the catchment. We will generate a stream network and catchment boundary using a digital terrain model and compare it with the actual stream network. Basic terrain statistics like slope will be calculated, along with a hypsometric curve to assess elevation distribution. Using the Thiessen polygon method, we will calculate areal precipitation and compare it to individual station data. Finally, land-use patterns will be analyzed using the 2012 Corine Land Cover (CLC) map to explore their impact on catchment hydrology.
# Import Catchment Area
catchment_area <- vect("basin/6210_Savinja_VSirje.shp")
catchment_area
## class : SpatVector
## geometry : polygons
## dimensions : 1, 2 (geometries, attributes)
## extent : 466905, 540775.8, 102666.7, 146214.3 (xmin, xmax, ymin, ymax)
## source : 6210_Savinja_VSirje.shp
## coord. ref. : gauss_krueger_SLO
## names : opomba AREA
## type : <chr> <num>
## values : NA 1.847e+09
# Import Digital Terrain Model (DTM)
digital_terrain_model <- rast("slo_100_ASCI.asc")
digital_terrain_model
## class : SpatRaster
## dimensions : 1625, 2478, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 375212, 623012, 30872, 193372 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source : slo_100_ASCI.asc
## name : slo_100_ASCI
# Plot Digital Terrain Model
plot(
digital_terrain_model,
main = "Digital Terrain Model (DEM)",
col = terrain.colors(20),
axes = TRUE
)
# Set the CRS for both the catchment area and the digital terrain model
crspost = crs(catchment_area)
crs(digital_terrain_model) = crspost
# Project the catchment area
project1 <- project(catchment_area, crspost)
# Plot the catchment area
plot(catchment_area, main = "Catchment Area Boundary", col = "lightblue", axes = TRUE)
# Overlay the projected catchment area
plot(project1, add = TRUE, border = "blue", lwd = 2)
# Add North Arrow and Scale Bar
north("topleft")
sbar(10000, xy = c(435000, 95000), type = "bar", divs = 2, below = "km", label = c(0, 5, 10))
# Load Actual Stream Network
actual_stream_network <- vect("streams/VT_POVR_CRTLine.shp")
# Fill Local Depressions
wbt_breach_depressions_least_cost(dem = "slo_100_ASCI.asc", output = "dem_breached.tif", dist = 500, fill = TRUE)
# Fill Depressions Using Wang and Liu Algorithm
wbt_fill_depressions_wang_and_liu(dem = "dem_breached.tif", output = "filled_dem.tif")
# Flow Accumulation Calculation
wbt_d8_flow_accumulation(input = "filled_dem.tif", output = "flow_accumulation.tif")
# Flow Direction Calculation
wbt_d8_pointer(dem = "filled_dem.tif", output = "flow_direction.tif")
# Generate River Network
wbt_extract_streams(flow_accum = "flow_accumulation.tif", output = "generated_streams.tif", threshold = 1000)
# Load Generated Stream Network
generated_stream_network <- rast("generated_streams.tif")
# Transform Generated Stream Network into Polygons
generated_stream_polygons <- as.polygons(generated_stream_network)
# Load Digital Terrain Model for Plotting
digital_terrain_model <- rast("slo_100_ASCI.asc")
# Compare Generated and Actual Stream Networks
plot(digital_terrain_model)
lines(generated_stream_polygons, col = "gray")
lines(actual_stream_network, col = "red")
legend("topleft", legend = c("Generated Streams", "Actual Streams"), col = c("gray", "red"),
lty = 1, lwd = 2)
plot(digital_terrain_model)
# Generate a catchment boundary
d1 = data.frame("veliko sirje", X=515395, Y=105435)
station = vect(d1, geom = c("X", "Y"))
points(station, lwd=100, col="red")
wbt_watershed(
d8_pntr = "flow_direction.tif",
pour_pts = "tocka.shp",
output = "catchment_area.tif"
)
# Import the catchment area
catchment_area_raster <- rast("catchment_area.tif")
# Transform raster to polygon
catchment_area_polygon <- as.polygons(catchment_area_raster)
# Plot the digital terrain model
plot(digital_terrain_model, main = "Catchment Area and Digital Terrain Model")
# Add the catchment area to the plot
plot(catchment_area_polygon, add = TRUE, col = "orange", border = "red")
lines(catchment_area, col = "white")
# Calculate Terrain Slope
terrain_slope <- terrain(digital_terrain_model, v = "slope", unit = "degrees")
# Calculate Basic Statistics
slope_statistics <- summary(values(terrain_slope))
slope_statistics
## slope
## Min. : 0.0
## 1st Qu.: 3.4
## Median : 8.9
## Mean :11.3
## 3rd Qu.:16.7
## Max. :72.0
## NA's :1888211
# Visualize Slope
plot(terrain_slope, col = terrain.colors(10), main = "Slope Map")
north("topleft")
# Histogram of Slope
hist(values(terrain_slope), main = "Slope Histogram", xlab = "Slope (degrees)", breaks = 20,
col = "lightblue")
The histogram represents the distribution of slope values within the Savinja Veliko Sirje catchment area. The majority of slopes are concentrated between 0° and 20°, with a significant peak in the lower range (0°–10°), indicating that the terrain is predominantly flat to gently sloping. Higher slope values (above 30°) are less frequent, suggesting that steep areas are limited in extent within the catchment. This slope distribution can have implications for runoff, erosion, and land-use planning within the region.
# Area and Perimeter of the Catchment
catchment_area_km2 <- expanse(catchment_area, unit = "km")
catchment_area_km2
## [1] 1845.736
catchment_perimeter <- perim(catchment_area)
catchment_perimeter
## [1] 272427.4
# Plot Hypsometric Curve
plot(ecdf(as.numeric(values(digital_terrain_model))), xlab = "Elevation (m)",
ylab = "Proportion of Area", main = "Hypsometric Curve of Catchment")
# Calculate and Display Elevation Quantiles
elevation_quantiles <- quantile(as.numeric(values(digital_terrain_model)), probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
abline(v = elevation_quantiles, col = c("red", "blue", "green"), lty = 2, lwd = 2)
The hypsometric curve shows the distribution of elevation across the catchment area. The majority of the catchment’s area is concentrated at lower elevations, as indicated by the steep initial rise in the curve. Elevation quantiles are marked for reference: 25% of the catchment lies below approximately 288 meters (red line), 50% lies below 466 meters (blue line), and 75% lies below 719 meters (green line). This distribution suggests a predominance of lower elevation areas, which can influence runoff dynamics, sediment transport, and vegetation distribution in the catchment.
elevation_quantiles
## 25% 50% 75%
## 287.6349 465.6508 719.4206
# whitebox::install_whitebox()
precip_data <- data.frame(names=c("P1_Luce", "P2_Solcava", "P3_Lasko", "P4_GornjiGrad", "P5_Radegunda"), lon=c(14.74, 14.69, 15.23, 14.8, 14.93), lat=c(46.35, 46.42, 46.15, 46.29, 46.36))
precip_data
## names lon lat
## 1 P1_Luce 14.74 46.35
## 2 P2_Solcava 14.69 46.42
## 3 P3_Lasko 15.23 46.15
## 4 P4_GornjiGrad 14.80 46.29
## 5 P5_Radegunda 14.93 46.36
slo_model <- elevation_30s(country="SI", path="slo_100_ASCI")
xprecip <- vect(precip_data, geom=c("lon", "lat"))
crs(xprecip) <- crs(slo_model)
plot(slo_model)
points(xprecip, col = "red", lwd = 10)
theisson <- voronoi(xprecip)
lines(theisson, col = "green")
# reproject catchment boundary to the lat-lon geographic coordinate system
catchment_area <- vect("basin/6210_Savinja_VSirje.shp")
catchment_theisson <- project(catchment_area, crs(slo_model))
lines(catchment_theisson, col = "blue")
# Intersection between the catchment area and thiessen polygons
theisson_area <- intersect(catchment_theisson, theisson)
lines(theisson_area, lwd=2, col="red")
# calculate the areas of polygons
expanse(theisson_area, unit = "km")
## [1] 160.6580 136.6530 923.7072 106.6205 518.0969
proportionP <- expanse(theisson_area, unit="km")/sum(expanse(theisson_area, unit="km"))
theisson_area$names
## [1] "P4_GornjiGrad" "P1_Luce" "P3_Lasko" "P2_Solcava"
## [5] "P5_Radegunda"
# areal precipitation
areal_precip <- main_data$P1_Luce*proportionP[2] + main_data$P2_Solcava*proportionP[4] + main_data$P3_Lasko*proportionP[3] + main_data$P4_GornjiGrad*proportionP[1] + main_data$P5_Radegunda*proportionP[5]
main_data$Areal_Precip <- areal_precip
plot(main_data$Date, areal_precip, type = "l", xlab = "Date", ylab = "Areal Precipitation [mm]")
# Load Land Cover map
land_use <- vect("land_use/CLC_2012_SIPolygon.shp")
# Plot the Land Cover map with the catchment area
plot(land_use, main = "Land Cover", axes = TRUE)
lines(catchment_area, col = "red")
# Check for invalid geometries in the land cover map and correct them
which(is.valid(land_use, TRUE) == FALSE)
## [1] 11139 14587 14588 14613 14674 14680 14689 14720 14770 14816 14829 14912
land_use <- makeValid(land_use)
# Perform spatial intersection
intersations <- intersect(catchment_area, land_use)
head(values(intersations))
## opomba AREA ID CODE_12 AREA_HA REMARK OPIS3
## 1 <NA> 1847142640 SI-106579 231 48.46294 <NA> Pašniki
## 2 <NA> 1847142640 SI-106609 231 48.02321 <NA> Pašniki
## 3 <NA> 1847142640 SI-106655 311 88.05743 <NA> Listnati gozd
## 4 <NA> 1847142640 SI-106427 231 131.10541 <NA> Pašniki
## 5 <NA> 1847142640 SI-106330 313 486.06507 <NA> Mešani gozd
## 6 <NA> 1847142640 SI-106435 231 70.65305 <NA> Pašniki
## LABEL3 OPIS2 LABEL2
## 1 Pastures Pašniki Pastures
## 2 Pastures Pašniki Pastures
## 3 Broad-leaved forest Gozdovi Forests
## 4 Pastures Pašniki Pastures
## 5 Mixed forest Gozdovi Forests
## 6 Pastures Pašniki Pastures
## OPIS1 LABEL1
## 1 Kmetijske povšine Agricultural areas
## 2 Kmetijske povšine Agricultural areas
## 3 Gozdne in deloma ohranjene naravne površine Forest and semi natural areas
## 4 Kmetijske povšine Agricultural areas
## 5 Gozdne in deloma ohranjene naravne površine Forest and semi natural areas
## 6 Kmetijske povšine Agricultural areas
## C_FAKTOR
## 1 0.004
## 2 0.004
## 3 0.002
## 4 0.004
## 5 0.002
## 6 0.004
# Plot the intersection
plot(intersations, "LABEL2", col = rainbow(20), plg = list(x = "bottomleft", cex = 0.7), mar = c(1, 1, 1, 1))
# Calculate the area of each land use
land_use_area = data.frame(area = expanse(intersations), desc = intersations$LABEL2)
# Calculate percentage for each land use type
land_use_area$percentage <- (land_use_area$area / sum(land_use_area$area)) * 100
head(land_use_area)
## area desc percentage
## 1 54342.355 Pastures 0.0029442111
## 2 54005.764 Pastures 0.0029259750
## 3 7562.658 Forests 0.0004097368
## 4 1147967.363 Pastures 0.0621956535
## 5 3512183.833 Forests 0.1902863929
## 6 669435.041 Pastures 0.0362692801
# Group land use by land use type
land_use_summary <- land_use_area %>%
group_by(desc) %>%
summarise(
total_area = sum(area),
percentage = (sum(area) / sum(land_use_area$area)) * 100
)
# Display summarized land use data
land_use_summary
## # A tibble: 12 × 3
## desc total_area percentage
## <chr> <dbl> <dbl>
## 1 Arable land 17402904. 0.943
## 2 Artificial, non-agricultural vegetated areas 2483004. 0.135
## 3 Forests 1032872691. 56.0
## 4 Heterogeneous agricultural areas 502325009. 27.2
## 5 Industrial, commercial and transport units 9315112. 0.505
## 6 Inland waters 5256309. 0.285
## 7 Mine, dump and construction sites 2329850. 0.126
## 8 Open spaces with little or no vegetation 20523663. 1.11
## 9 Pastures 161496582. 8.75
## 10 Permanent crops 17182187. 0.931
## 11 Scrub and/or herbaceous vegetation associations 39678600. 2.15
## 12 Urban fabric 34869761. 1.89
land_use_colors <- rainbow(nrow(land_use_summary))
# Prepare data for plotting
land_use_summary$desc <- factor(land_use_summary$desc, levels = land_use_summary$desc)
# Create the plot
ggplot(land_use_summary, aes(x = percentage, y = desc, fill = desc)) +
geom_bar(stat = "identity", color = "black", width = 0.5) +
scale_fill_manual(values = land_use_colors) +
labs(
title = "Land Use Percentages in Catchment Area",
x = "Percentage (%)",
y = ""
) +
theme_minimal() +
theme(
legend.position = "none",
axis.text.y = element_text(size = 6, hjust = 1, color = "black"),
axis.text.x = element_text(size = 8),
plot.title = element_text(size = 14, face = "bold")
) +
geom_text(
aes(label = paste0(round(percentage, 2), "%")),
hjust = -0.2,
color = "black",
size = 3
) +
coord_cartesian(clip = "off")
The bar chart represents the distribution of land use types within the catchment area. Forests dominate the catchment, covering 55.96% of the total area, followed by heterogeneous agricultural areas at 27.22%. Pastures account for 8.75%, while urban fabric and scrub or herbaceous vegetation associations represent 1.89% and 2.15%, respectively. Other land use types, including permanent crops, inland waters, and industrial or transport units, each occupy less than 1% of the catchment.
In this section, we prepare, calibrate, validate, and compare hydrological models for the catchment. The GR4J and CemaNeige GR6J models will be used to simulate catchment hydrology. Relevant data will be prepared for both models, and datasets will be split (50% for calibration and validation) with the first year reserved as a warm-up period. Performance metrics of the models will be evaluated and compared to determine their suitability for modeling the catchment’s hydrology.
# Create data frame
model_data <- main_data[, c("Date", "Discharge", "Runoff", "Areal_Precip", "Evapotrans", "Airtemp")]
# Display the new dataframe
head(model_data)
## Date Discharge Runoff Areal_Precip Evapotrans Airtemp
## 1 2010-01-01 59.058 2.762648 3.9654273 0.3 5.7
## 2 2010-01-02 53.925 2.522534 2.1476618 0.9 3.5
## 3 2010-01-03 47.681 2.230449 0.9421236 0.5 -2.0
## 4 2010-01-04 41.040 1.919792 0.0000000 0.2 -3.0
## 5 2010-01-05 38.694 1.810050 5.2556681 0.2 -1.8
## 6 2010-01-06 36.966 1.729216 7.3304657 0.2 -1.0
# Display basic statistics
summary(model_data)
## Date Discharge Runoff
## Min. :2010-01-01 00:00:00.00 Min. : 5.515 Min. : 0.2580
## 1st Qu.:2013-04-01 18:00:00.00 1st Qu.: 13.966 1st Qu.: 0.6533
## Median :2016-07-01 12:00:00.00 Median : 22.032 Median : 1.0306
## Mean :2016-07-01 12:24:47.61 Mean : 38.480 Mean : 1.8000
## 3rd Qu.:2019-10-01 06:00:00.00 3rd Qu.: 40.252 3rd Qu.: 1.8829
## Max. :2022-12-31 00:00:00.00 Max. :859.744 Max. :40.2176
## Areal_Precip Evapotrans Airtemp
## Min. : 0.00000 Min. :0.000 Min. :-13.80
## 1st Qu.: 0.00000 1st Qu.:0.700 1st Qu.: 4.40
## Median : 0.01672 Median :1.800 Median : 11.30
## Mean : 3.60893 Mean :2.264 Mean : 10.88
## 3rd Qu.: 2.57378 3rd Qu.:3.600 3rd Qu.: 17.70
## Max. :109.52009 Max. :7.700 Max. : 28.90
# Plot flow data
plot(model_data$Runoff, type = "l", col = "blue", lwd = 2, main = "Flow Data Over Time", xlab = "Days", ylab = "Flow (mm)")
# Create Inputs Model
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = as.POSIXct(model_data$Date, format = "%Y-%m-%d"), Precip = model_data$Areal_Precip, PotEvap = model_data$Evapotrans)
# Identify the index where calibration ends
calibration_end <- which(format(as.Date(model_data$Date), format = "%Y-%m-%d") == "2016-12-31")
# Define the sequence for the run period
Ind_Run <- seq(366, calibration_end)
# Create RunOptions with defined run and warm-up periods
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run, IndPeriod_WarmUp = 1:365)
# Extract observed discharge data for the run period
Obs_Discharge <- model_data$Runoff[Ind_Run]
# Create InputsCrit using the Nash-Sutcliffe Efficiency (NSE) criterion
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, RunOptions = RunOptions, Obs = Obs_Discharge)
# Define calibration options using the Michel algorithm
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
# Perform calibration to find optimal model parameters
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit, CalibOptions = CalibOptions, FUN_MOD = RunModel_GR4J)
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
## Screening completed (81 runs)
## Param = 432.681, -0.020, 42.098, 1.417
## Crit. NSE[Q] = 0.5850
## Steepest-descent local search in progress
## Calibration completed (27 iterations, 282 runs)
## Param = 549.741, 0.270, 45.955, 0.548
## Crit. NSE[Q] = 0.8192
# Extract calibrated parameter values
ParamGR4J <- OutputsCalib$ParamFinalR
ParamGR4J
## [1] 549.7413151 0.2699136 45.9554779 0.5477049
# Run the GR4J model with calibrated parameters
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = ParamGR4J)
# Plot the comparison between measured and modeled data
plot(OutputsModel, model_data$Runoff[366:calibration_end])
# Define the sequence for the validation period
Ind_Run1 <- seq((calibration_end + 1), nrow(model_data))
# Create RunOptions for the validation period
RunOptions1 <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, IndPeriod_Run = Ind_Run1)
## Warning in CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, : model warm up period not defined: default configuration used
## the year preceding the run period is used
# Run the GR4J model for the validation period
OutputsModel1 <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions1, Param = ParamGR4J)
# Plot the comparison between measured and modeled runoff for validation
plot(OutputsModel1, model_data$Runoff[Ind_Run1])
# Create InputsCrit for validation using NSE as the criterion
InputsCrit1 <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, RunOptions = RunOptions1, Obs = model_data$Runoff[Ind_Run1])
# Calculate the validation performance (NSE)
OutputsCrit1 <- ErrorCrit_NSE(InputsCrit = InputsCrit1, OutputsModel = OutputsModel1)
## Crit. NSE[Q] = 0.7192
# Print the validation results
OutputsCrit1
## $CritValue
## [1] 0.7191667
##
## $CritName
## [1] "NSE[Q]"
##
## $CritBestValue
## [1] 1
##
## $Multiplier
## [1] -1
##
## $Ind_notcomputed
## NULL
##
## attr(,"class")
## [1] "NSE" "ErrorCrit"
# Prepare observation data in the required format
preob <- data.frame(DatesR = as.POSIXct(model_data$Date, format = "%Y-%m-%d"), P=model_data[,4], E=model_data[, 5],Qmm=model_data[,3])
# Display the first few rows of the observation data
head(preob)
## DatesR P E Qmm
## 1 2010-01-01 3.9654273 0.3 2.762648
## 2 2010-01-02 2.1476618 0.9 2.522534
## 3 2010-01-03 0.9421236 0.5 2.230449
## 4 2010-01-04 0.0000000 0.2 1.919792
## 5 2010-01-05 5.2556681 0.2 1.810050
## 6 2010-01-06 7.3304657 0.2 1.729216
# Load example data for the basin
data(L0123001)
# Display the hypsometric curve data format
BasinInfo
## $BasinCode
## [1] "L0123001"
##
## $BasinName
## [1] "Blue River at Nourlangie Rock"
##
## $BasinArea
## [1] 360
##
## $HypsoData
## [1] 286 309 320 327 333 338 342 347 351 356 360 365 369 373 378
## [16] 382 387 393 399 405 411 417 423 428 434 439 443 448 453 458
## [31] 463 469 474 480 485 491 496 501 507 513 519 524 530 536 542
## [46] 548 554 560 566 571 577 583 590 596 603 609 615 622 629 636
## [61] 642 649 656 663 669 677 684 691 698 706 714 722 730 738 746
## [76] 754 762 770 777 786 797 808 819 829 841 852 863 875 887 901
## [91] 916 934 952 972 994 1012 1029 1054 1080 1125 1278
# Calculate the hypsometric curve using the digital terrain model
hypso <- c(min(as.numeric(values(digital_terrain_model)),na.rm=T),as.numeric(quantile(as.numeric(values(digital_terrain_model)),probs=seq(from=0.01,by=0.01,to=0.99),na.rm=T)),max(as.numeric(values(digital_terrain_model)),na.rm=T))
# Define model settings, including the number of height zones and the hypsometric curve
InputsModel2 <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR6J, as.POSIXct(model_data$Date, format = "%Y-%m-%d"), Precip = model_data[,4], PotEvap = model_data[, 5], TempMean = model_data[,6], ZInputs = median(hypso), HypsoData = hypso, NLayers = 5)
## input series were successfully created on 5 elevation layers for use by CemaNeige
# Define run options, including warm-up period and annual solid precipitation
RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, IndPeriod_Run = Ind_Run, IndPeriod_WarmUp = 1:365, MeanAnSolidPrecip=c(70,130,180,250,350))
# Define the goodness of fit criterion and discharge data
InputsCrit2 <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel2, RunOptions = RunOptions2,
Obs = model_data[366:calibration_end,3])
# Define the calibration method
CalibOptions2 <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR6J, FUN_CALIB = Calibration_Michel)
# Run the model calibration
OutputsCalib2 <- Calibration_Michel(InputsModel = InputsModel2, RunOptions = RunOptions2, InputsCrit = InputsCrit2, CalibOptions = CalibOptions2, FUN_MOD = RunModel_CemaNeigeGR6J)
## Grid-Screening in progress (
## 0% 20% 40% 60% 80% 100%)
## Screening completed (6561 runs)
## Param = 90.017, -0.521, 27.113, 1.369, 0.220, 54.598, 0.002, 6.764
## Crit. NSE[Q] = 0.6486
## Steepest-descent local search in progress
## Calibration completed (65 iterations, 7565 runs)
## Param = 178.520, 0.107, 14.540, 1.042, 0.447, 32.407, 0.002, 6.540
## Crit. NSE[Q] = 0.8876
# Save the calibrated parameter values
ParamGR6JCemaNeige <- OutputsCalib2$ParamFinalR
# Run the model with the calibrated parameter values
OutputsModel2 <- RunModel_CemaNeigeGR6J(InputsModel = InputsModel2, RunOptions = RunOptions2,
Param = ParamGR6JCemaNeige)
# Plot the observed versus simulated values
plot(OutputsModel2, model_data[366:calibration_end,3])
# Calculate NSE for the validation period
ErrorCrit_NSE(InputsCrit = InputsCrit2, OutputsModel = OutputsModel2)
## Crit. NSE[Q] = 0.8876
## $CritValue
## [1] 0.8876181
##
## $CritName
## [1] "NSE[Q]"
##
## $CritBestValue
## [1] 1
##
## $Multiplier
## [1] -1
##
## $Ind_notcomputed
## NULL
##
## attr(,"class")
## [1] "NSE" "ErrorCrit"
# Create run options for validation
RunOptions3 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, IndPeriod_Run = Ind_Run1)
## Warning in CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, : model warm up period not defined: default configuration used
## the year preceding the run period is used
## Warning in CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel2, : 'MeanAnSolidPrecip' not defined: it was automatically set to c(146,146,146,146,146) from the 'InputsModel' given to the function. Be careful in case your application is (short-term) forecasting.
# Run the model for validation
OutputsModel3 <- RunModel_CemaNeigeGR6J(InputsModel = InputsModel2, RunOptions = RunOptions3, Param = ParamGR6JCemaNeige)
# Plot simulated vs observed discharge
plot(OutputsModel3, Qobs = model_data[Ind_Run1,3])
# Create inputs for the NSE criterion during validation
InputsCrit3 <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel2, RunOptions = RunOptions3, Obs = model_data[Ind_Run1,3])
# Calculate NSE for the validation period
OutputsCrit3 <- ErrorCrit_NSE(InputsCrit = InputsCrit3, OutputsModel = OutputsModel3)
## Crit. NSE[Q] = 0.7892
# Calculate NSE for GR4J model
OutputsCrit1 <- ErrorCrit_NSE(InputsCrit = InputsCrit1, OutputsModel = OutputsModel1)
## Crit. NSE[Q] = 0.7192
NSE_GR4J <- OutputsCrit1$CritValue
# Calculate NSE for GR6J model
OutputsCrit3 <- ErrorCrit_NSE(InputsCrit = InputsCrit3, OutputsModel = OutputsModel3)
## Crit. NSE[Q] = 0.7892
NSE_GR6J <- OutputsCrit3$CritValue
# Compare the performance of the two models
if (NSE_GR6J > NSE_GR4J) {
cat("The GR6J model performs better with NSE =", NSE_GR6J,
"compared to GR4J with NSE =", NSE_GR4J, "\n")
} else if (NSE_GR6J < NSE_GR4J) {
cat("The GR4J model performs better with NSE =", NSE_GR4J,
"compared to GR6J with NSE =", NSE_GR6J, "\n")
} else {
cat("Both models perform equally well with NSE =", NSE_GR6J, "\n")
}
## The GR6J model performs better with NSE = 0.7891693 compared to GR4J with NSE = 0.7191667
The comparison of model performance based on Nash-Sutcliffe Efficiency (NSE) indicates that the GR6J model outperforms the GR4J model. The GR6J model achieved an NSE value of 0.789, whereas the GR4J model achieved an NSE value of 0.719. The enhanced performance of the GR6J model can be attributed to its additional parameters, including an exponential store that improves low-flow simulations, allowing it to better capture hydrological processes compared to the GR4J model. This suggests that GR6J provides a more accurate representation of observed streamflow dynamics.
In this section, we will explore hydrological modeling using the TUWIEN model. The analysis begins by preparing all relevant data for the model. Following data preparation, the TUWIEN model will be calibrated and validated by splitting the dataset into training and validation periods, using the first year of data for warm-up. Finally, the performance of the TUWIEN model will be compared with the CemaNeige GR6J model. Both graphical and statistical metrics will be used to assess the models’ capabilities.
# Create a data frame for the TUWIEN model
preob1 <- data.frame(
DatesR = as.POSIXct(model_data$Date, format = "%Y-%m-%d"),
P = model_data[, 4],
E = model_data[, 5],
Qmm = model_data[, 3],
T = model_data[, 6]
)
head(preob1)
## DatesR P E Qmm T
## 1 2010-01-01 3.9654273 0.3 2.762648 5.7
## 2 2010-01-02 2.1476618 0.9 2.522534 3.5
## 3 2010-01-03 0.9421236 0.5 2.230449 -2.0
## 4 2010-01-04 0.0000000 0.2 1.919792 -3.0
## 5 2010-01-05 5.2556681 0.2 1.810050 -1.8
## 6 2010-01-06 7.3304657 0.2 1.729216 -1.0
# Plot precipitation and runoff data
plot(preob1$DatesR, preob1$P, type = "l", col = "blue", ylab = "Precipitation [mm]", xlab = "Date", main = "Precipitation over Time")
plot(preob1$DatesR, preob1$Qmm, type = "l", col = "green", ylab = "Runoff [mm]", xlab = "Date", main = "Runoff over Time")
# End of calibration period
konec <- which(format(preob1$DatesR, format = "%Y-%m-%d") == "2016-12-31")
# Run TUWIEN model with randomly selected parameter values
sim1 <- TUWmodel(
prec = preob1$P[1:konec],
airt = preob1$T[1:konec],
ep = preob1$E[1:konec],
area = 1,
param = c(1.3, 2.0, -1.0, 1.0, 0.0, 0.8, 360.0, 0.2, 0.3, 7.0, 150.0, 50.0, 2.0, 10.0, 25.0),
incon = c(50, 0, 2.5, 2.5)
)
# Check the structure of the simulation output
str(sim1)
## List of 22
## $ itsteps: int 2557
## $ nzones : int 1
## $ area : num 1
## $ param : num [1, 1:15] 1.3 2 -1 1 0 0.8 360 0.2 0.3 7 ...
## ..- attr(*, "names")= chr [1:15] "SCF" "DDF" "Tr" "Ts" ...
## $ incon : num [1, 1:4] 50 0 2.5 2.5
## ..- attr(*, "names")= chr [1:4] "SSM0" "SWE0" "SUZ0" "SLZ0"
## $ prec : num [1:2557] 3.965 2.148 0.942 0 5.256 ...
## $ airt : num [1:2557] 5.7 3.5 -2 -3 -1.8 -1 -0.3 -0.4 0.6 0.9 ...
## $ ep : num [1:2557] 0.3 0.9 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...
## $ output : num [1, 1:20, 1:2557] 0.404 0 51.24 3.965 0 ...
## $ qzones : num [1, 1:2557] 0.40409 0.1516 0.15334 0.00696 0.01563 ...
## $ q : num [1, 1:2557] 0.40409 0.1516 0.15334 0.00696 0.01563 ...
## $ swe : num [1, 1:2557] 0 0 1.22 1.22 8.06 ...
## $ melt : num [1, 1:2557] 0 0 0 0 0 0 0 0 1.2 1.8 ...
## $ q0 : num [1, 1:2557] 0 0 0 0 0 0 0 0 0 0 ...
## $ q1 : num [1, 1:2557] 0.374 0.26 0 0 0 ...
## $ q2 : num [1, 1:2557] 0.0298 0.0429 0.0558 0.0554 0.055 ...
## $ moist : num [1, 1:2557] 51.2 51.8 51.8 51.8 51.8 ...
## $ rain : num [1, 1:2557] 3.97 2.15 0 0 0 ...
## $ snow : num [1, 1:2557] 0 0 0.942 0 5.256 ...
## $ eta : num [1, 1:2557] 0.0534 0.1623 0 0 0 ...
## $ suz : num [1, 1:2557] 2.8 1.99 0 0 0 ...
## $ slz : num [1, 1:2557] 4.47 6.43 8.36 8.31 8.25 ...
# Plot simulated and observed runoff for the calibration period
plot(as.numeric(sim1$q)[366:konec], type = "l", col = "red", ylab = "Discharge [mm]", xlab = "Days", main = "Calibration Results")
lines(preob1$Qmm[366:konec], col = "green")
legend("topright", legend = c("Simulated", "Observed"), col = c("red", "green"), lty = 1)
# plot the simulated values (for the first year)
plot(as.numeric(sim1$q)[1:365],type="l", col="red",ylab="Discharge [mm]")
# add measured flow data
lines(preob1[1:365,4], col="green")
# Calculate goodness-of-fit statistics for the calibration period
gof_results <- gof(as.numeric(sim1$q)[366:konec], preob1$Qmm[366:konec])
gof_results
## [,1]
## ME 1.11
## MAE 1.45
## MSE 3.67
## RMSE 1.92
## ubRMSE 1.56
## NRMSE % 76.30
## PBIAS % 59.80
## RSR 0.76
## rSD 0.86
## NSE 0.42
## mNSE -0.02
## rNSE -0.74
## wNSE 0.74
## wsNSE 0.27
## d 0.83
## dr 0.49
## md 0.50
## rd 0.49
## cp -0.09
## r 0.79
## R2 0.42
## bR2 0.42
## VE 0.22
## KGE 0.35
## KGElf 0.20
## KGEnp 0.33
## KGEkm 0.32
# Define the objective function for optimization
msespr <- function (param, precip, temp, potevap, runoff, area) {
simu <- TUWmodel(param, prec=as.numeric(precip), airt=as.numeric(temp), ep=as.numeric(potevap), area=area)$q
simu <- simu[-c(1:365)] # Remove warm-up period
obse <- runoff[-c(1:365)]
mse(simu, obse) # Calculate mean squared error
}
# Calibrate the model using DEoptim
calibrate_period1 <- DEoptim(
fn = msespr,
lower = c(0.9, 0.0, 1.0, -3.0, -2.0, 0.0, 0.0, 0.0, 0.0, 2.0, 30.0, 1.0, 0.0, 0.0, 0.0),
upper = c(1.5, 5.0, 3.0, 1.0, 2.0, 1.0, 600.0, 20.0, 2.0, 30.0, 250.0, 100.0, 8.0, 30.0, 50.0),
control = DEoptim.control(itermax = 60),
precip = preob1$P[1:konec],
temp = preob1$T[1:konec],
potevap = preob1$E[1:konec],
runoff = preob1$Qmm[1:konec],
area = 1
)
## Iteration: 1 bestvalit: 1.562010 bestmemit: 1.258932 0.834634 2.819617 -0.828915 -1.962681 0.986918 319.547874 17.270585 1.128826 7.168388 177.013609 22.856899 3.719463 2.875780 40.803773
## Iteration: 2 bestvalit: 1.562010 bestmemit: 1.258932 0.834634 2.819617 -0.828915 -1.962681 0.986918 319.547874 17.270585 1.128826 7.168388 177.013609 22.856899 3.719463 9.268427 40.803773
## Iteration: 3 bestvalit: 1.515548 bestmemit: 1.134806 4.458489 2.786875 -2.867137 0.428393 0.875656 240.440627 19.381912 0.676021 26.449547 107.558175 3.290968 6.649246 23.730181 17.794380
## Iteration: 4 bestvalit: 1.453108 bestmemit: 1.053977 2.367076 1.982698 -1.581068 -1.023994 0.672748 105.929010 17.989298 1.721246 7.707436 168.114704 29.029833 3.775845 16.283480 11.133193
## Iteration: 5 bestvalit: 1.264874 bestmemit: 1.105428 0.707724 1.829851 0.186863 -0.934620 0.946221 279.553487 8.838589 1.762308 8.537585 38.952472 4.706629 5.459383 23.533203 46.675108
## Iteration: 6 bestvalit: 1.264874 bestmemit: 1.105428 0.707724 1.829851 0.186863 -0.934620 0.946221 279.553487 8.838589 1.762308 8.537585 38.952472 4.706629 5.459383 23.533203 46.675108
## Iteration: 7 bestvalit: 1.236245 bestmemit: 1.430549 0.786743 2.201725 -2.233606 -0.381786 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 10.711176 49.275859
## Iteration: 8 bestvalit: 1.231692 bestmemit: 1.430549 0.786743 2.006815 -2.233606 -0.381786 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 10.711176 49.275859
## Iteration: 9 bestvalit: 1.231692 bestmemit: 1.430549 0.786743 2.006815 -2.233606 -0.381786 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 10.711176 49.275859
## Iteration: 10 bestvalit: 1.198061 bestmemit: 0.943814 1.478951 1.403058 -1.183636 -1.956553 0.700740 315.402509 1.717766 1.238394 3.304905 151.668674 23.085526 1.803697 28.368848 49.664354
## Iteration: 11 bestvalit: 1.198061 bestmemit: 0.943814 1.478951 1.403058 -1.183636 -1.956553 0.700740 315.402509 1.717766 1.238394 3.304905 151.668674 23.085526 1.803697 28.368848 49.664354
## Iteration: 12 bestvalit: 1.198061 bestmemit: 0.943814 1.478951 1.403058 -1.183636 -1.956553 0.700740 315.402509 1.717766 1.238394 3.304905 151.668674 23.085526 1.803697 28.368848 49.664354
## Iteration: 13 bestvalit: 1.198061 bestmemit: 0.943814 1.478951 1.403058 -1.183636 -1.956553 0.700740 315.402509 1.717766 1.238394 3.304905 151.668674 23.085526 1.803697 28.368848 49.664354
## Iteration: 14 bestvalit: 1.198061 bestmemit: 0.943814 1.478951 1.403058 -1.183636 -1.956553 0.700740 315.402509 1.717766 1.238394 3.304905 151.668674 23.085526 1.803697 28.368848 49.664354
## Iteration: 15 bestvalit: 0.985166 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.381786 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 10.711176 49.275859
## Iteration: 16 bestvalit: 0.985166 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.381786 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 17 bestvalit: 0.985166 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.381786 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 18 bestvalit: 0.985061 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 19 bestvalit: 0.985061 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 20 bestvalit: 0.985061 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 21 bestvalit: 0.985061 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 22 bestvalit: 0.985061 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 14.606405 2.894247 0.479086 32.196982
## Iteration: 23 bestvalit: 0.965310 bestmemit: 0.948157 1.080407 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 24 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 25 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 26 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 27 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 28 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 29 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 30 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 31 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 32 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 33 bestvalit: 0.949514 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.813976 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 34 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 35 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 36 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 37 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 38 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 39 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 40 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 41 bestvalit: 0.946884 bestmemit: 0.948157 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 42 bestvalit: 0.936877 bestmemit: 0.922315 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 43 bestvalit: 0.936877 bestmemit: 0.922315 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 44 bestvalit: 0.936877 bestmemit: 0.922315 1.392848 2.006815 -2.233606 -0.269818 0.998964 352.814081 2.247926 1.932069 2.156076 40.920631 16.587395 2.894247 0.479086 32.196982
## Iteration: 45 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 46 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 47 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 48 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 49 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 50 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 51 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 52 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 53 bestvalit: 0.935961 bestmemit: 0.960670 1.386594 2.115617 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 54 bestvalit: 0.935121 bestmemit: 0.960670 1.386594 1.892064 0.171584 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 55 bestvalit: 0.933865 bestmemit: 0.960670 1.386594 1.892064 0.330761 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 56 bestvalit: 0.933865 bestmemit: 0.960670 1.386594 1.892064 0.330761 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 6.521182 4.303730 8.403318 33.161193
## Iteration: 57 bestvalit: 0.895191 bestmemit: 0.960670 1.386594 1.892064 0.330761 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 3.314976 4.303730 8.403318 33.161193
## Iteration: 58 bestvalit: 0.895191 bestmemit: 0.960670 1.386594 1.892064 0.330761 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 3.314976 4.303730 8.403318 33.161193
## Iteration: 59 bestvalit: 0.895191 bestmemit: 0.960670 1.386594 1.892064 0.330761 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 3.314976 4.303730 8.403318 33.161193
## Iteration: 60 bestvalit: 0.895191 bestmemit: 0.960670 1.386594 1.892064 0.330761 -1.342859 0.958450 297.556700 3.363709 1.763296 7.986158 34.336170 3.314976 4.303730 8.403318 33.161193
# Run the model with calibrated parameters
sim2 <- TUWmodel(
prec = preob1$P[1:konec],
airt = preob1$T[1:konec],
ep = preob1$E[1:konec],
area = 1,
param = calibrate_period1$optim$bestmem
)
# Plot simulated vs observed runoff for the validation period
plot(as.numeric(sim2$q)[366:konec], type = "l", col = "red", ylab = "Discharge [mm]", xlab = "Days", main = "Validation Results")
lines(preob1$Qmm[366:konec], col = "green")
legend("topright", legend = c("Simulated", "Observed"), col = c("red", "green"), lty = 1)
# Plot simulated vs observed runoff for the first year
plot(as.numeric(sim2$q)[1:365], type = "l", col = "red", ylab = "Discharge [mm]", xlab = "Days", main = "Validation Results for First Year")
lines(preob1[1:365, 4], col = "green")
# Calculate NSE for the validation period
NSE_validation <- NSE(as.numeric(sim2$q)[366:konec], preob1$Qmm[366:konec])
print(paste("NSE for validation period: ", round(NSE_validation, 3)))
## [1] "NSE for validation period: 0.858"
# Calculate NSE for GR6J model
OutputsCrit3 <- ErrorCrit_NSE(InputsCrit = InputsCrit3, OutputsModel = OutputsModel3)
## Crit. NSE[Q] = 0.7892
NSE_GR6J <- OutputsCrit3$CritValue
# Calculate NSE for TUWmodel
NSE_TUW <- NSE(as.numeric(sim2$q)[366:konec], preob1$Qmm[366:konec])
# Compare the performance of the two models
if (NSE_GR6J > NSE_TUW) {
cat("The GR6J model performs better with NSE =", NSE_GR6J,
"compared to TUWmodel with NSE =", NSE_TUW, "\n")
} else if (NSE_GR6J < NSE_TUW) {
cat("The TUWmodel performs better with NSE =", NSE_TUW,
"compared to GR6J with NSE =", NSE_GR6J, "\n")
} else {
cat("Both models perform equally well with NSE =", NSE_GR6J, "\n")
}
## The TUWmodel performs better with NSE = 0.8579348 compared to GR6J with NSE = 0.7891693
The comparison of model performance based on Nash-Sutcliffe Efficiency (NSE) reveals that the TUWmodel outperforms the GR6J model. The TUWmodel achieved an NSE value of 0.858, compared to the GR6J model’s NSE value of 0.789. The superior performance of the TUWmodel can be attributed to its more sophisticated representation of hydrological processes, potentially including enhanced snowmelt modeling and soil moisture dynamics. This indicates that the TUWmodel provides a more precise simulation of observed streamflow, particularly in regions or conditions where such processes are significant.
In this section, we explore the use of stochastic simulation techniques to model precipitation. First, a stochastic precipitation simulator will be fitted to measured data using the GWEX package. Ten realizations of 13-year simulations will then be generated and compared to the measured data to evaluate the simulator’s performance. Finally, the simulated precipitation data will serve as input for the previously calibrated hydrological model (CemaNeige GR6J) to assess its performance under stochastic conditions.
# Transform Precipitation Data to Matrix
precip_matrix <- as.matrix(main_data[, 8:12])
# Define Precipitation Observations for GWEX
observed_precip <- GwexObs(variable = 'Prec', date = as.Date(main_data[, 1], format = "%Y.%m.%d"),
obs = precip_matrix)
# Estimate Parameters for Precipitation Simulator
precip_params <- fitGwexModel(observed_precip, listOption = list(th = 0.5, nChainFit = 1000))
## [1] "Fit generator"
precip_params
## *** Class GWex ***
## * type of variable = Prec
## * Version = v1.1.0
##
## #### Options for the occurrence process ####
## threshold: 0.5
## order of the Markov chain: 2
##
## #### Options for the amount process ####
## spatial dependence: Gaussian copula
## apply a MAR(1) process: FALSE
## apply a master Markov chain:
# Simulate Precipitation (10 Realizations for 13 Years)
simulated_precip <- simGwexModel(precip_params, nb.rep = 10, d.start = as.Date("01012025", "%d%m%Y"),
d.end = as.Date("31122037", "%d%m%Y"))
## [1] "Generate scenarios"
# Summary of Simulated Precipitation
precip_summary <- summary(simulated_precip@sim[, 1, 1])
precip_summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 4.089 2.486 183.289
P1_Luce_summary <- summary(main_data$P1_Luce)
P1_Luce_summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 4.236 2.200 116.100
# Plot Simulated Precipitation for Station 1, Realization 1
plot(simulated_precip@date, simulated_precip@sim[,1,1], type = "l", col = "green",
main = "Simulated Precipitation Station 1 (P1_Luce)", xlab = "Date",
ylab = "Precipitation (mm)")
# Add Realization 2
lines(simulated_precip@date, simulated_precip@sim[,1,2], col = "purple")
legend("topright", legend = c("Realization 1", "Realization 2"), col = c("green", "purple"), lty = 1)
# P1_Luce measured yearly precipitation
P1_Luce_yearly <- main_data %>%
group_by(Year) %>%
summarise(P1_Luce_measured = sum(P1_Luce))
P1_Luce_yearly
## # A tibble: 13 × 2
## Year P1_Luce_measured
## <dbl> <dbl>
## 1 2010 1726.
## 2 2011 1338.
## 3 2012 1712.
## 4 2013 1538.
## 5 2014 2088.
## 6 2015 1261
## 7 2016 1610.
## 8 2017 1652.
## 9 2018 1436.
## 10 2019 1553
## 11 2020 1503.
## 12 2021 1554.
## 13 2022 1139.
# Dataframe of dates and years
simulated_yearly <- data.frame(Date = simulated_precip@date,
Year = format(simulated_precip@date, "%Y"))
# Add yearly precipitation for each realization
for (i in 1:10) {
simulated_yearly[[paste0("Realization_", i)]] <- simulated_precip@sim[, 1, i]
}
# Calculate yearly totals for each realization
simulated_combined <- simulated_yearly %>%
group_by(Year) %>%
summarise(across(starts_with("Realization"), ~ sum(.x, na.rm = TRUE)))
simulated_combined
## # A tibble: 13 × 11
## Year Realization_1 Realization_2 Realization_3 Realization_4 Realization_5
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2025 1207. 1623. 1736. 1846. 1483.
## 2 2026 1541. 1417. 1396. 1534. 1476.
## 3 2027 1810. 1613. 1650. 1752. 1652.
## 4 2028 1431. 1558. 1804. 1185. 1772.
## 5 2029 1256. 1774. 1728. 1343. 1429.
## 6 2030 1401. 1410. 1427. 1524. 1562.
## 7 2031 1495. 1792. 1685. 1455. 1734.
## 8 2032 1590. 1108. 1183. 1577. 1627.
## 9 2033 1226. 1349. 1787. 2001. 1268.
## 10 2034 1658. 1241. 1922. 1388. 1082.
## 11 2035 1275. 1258. 1308. 1480. 1672.
## 12 2036 1964. 1550. 1705. 1535. 1620.
## 13 2037 1562. 1613. 1198. 1757. 1291.
## # ℹ 5 more variables: Realization_6 <dbl>, Realization_7 <dbl>,
## # Realization_8 <dbl>, Realization_9 <dbl>, Realization_10 <dbl>
# Convert Year to numeric
P1_Luce_yearly <- P1_Luce_yearly %>% mutate(Year = as.numeric(Year))
simulated_combined <- simulated_combined %>% mutate(Year = as.numeric(Year))
# Set up the plot for measured precipitation
plot(P1_Luce_yearly$Year, P1_Luce_yearly$P1_Luce_measured, type = "o", col = "blue", pch = 16, lwd = 2,
xlab = "Year", ylab = "Total Precipitation (mm)", main = "Measured vs Simulated Yearly Precipitation(P1_Luce)",
ylim = c(500, 2500))
# Add simulated precipitation for each realization
for (i in 1:10) {
lines(P1_Luce_yearly$Year, simulated_combined[[paste0("Realization_", i)]],
col = "red", lty = 2)
}
# Add a legend
legend("topright", legend = c("Measured", "Simulated (Realizations)"),
col = c("blue", "red"), lty = c(1, 2), pch = c(16, NA), lwd = c(2, 1))
for (i in 1:10) {
calibration_end1 <- which(format(as.Date(simulated_precip@date), format = "%Y-%m-%d") == "2031-12-31")
Ind_Run2 <- seq(366, calibration_end1)
InputsModel3 <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR6J,
DatesR = as.POSIXct(simulated_precip@date, format = "%Y-%m-%d"),
Precip = simulated_precip@sim[, 1, i], PotEvap = model_data[, 5],
TempMean = model_data[, 6], ZInputs = median(hypso),
HypsoData = hypso, NLayers = 5)
RunOptions3 <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel3,
IndPeriod_Run = Ind_Run2, IndPeriod_WarmUp = 1:365,
MeanAnSolidPrecip = c(70, 130, 180, 250, 350))
OutputsModel3 <- RunModel_CemaNeigeGR6J(InputsModel = InputsModel3, RunOptions = RunOptions3,
Param = ParamGR6JCemaNeige)
plot(OutputsModel3)
}
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
## input series were successfully created on 5 elevation layers for use by CemaNeige
## Warning in plot.OutputsModel(OutputsModel3): the "Error" and "CorQQ" plot(s)
## cannot be drawn if there is no 'Qobs'
In this section, we analyze the surface air temperature data from MODIS for the Savinja_VSirje catchment over a period of at least two years. The temperature characteristics are examined and compared with measured discharge data to identify potential correlations and impacts on hydrological behavior. Additionally, ERA5 data for the entire available period is downloaded and analyzed. A comparison between ERA5 spatial precipitation data and gauge-based areal precipitation measurements is conducted, with plots provided to visually compare the datasets and highlight any differences or similarities.
# List all available MODIS products
head(mt_products())
## product
## 1 Daymet
## 2 ECO4ESIPTJPL
## 3 ECO4WUE
## 4 GEDI03
## 5 GEDI04_B
## 6 MCD12Q1
## description
## 1 Daily Surface Weather Data (Daymet) on a 1-km Grid for North America, Version 4 R1
## 2 ECOSTRESS Evaporative Stress Index PT-JPL (ESI) Daily L4 Global 70 m
## 3 ECOSTRESS Water Use Efficiency (WUE) Daily L4 Global 70 m
## 4 GEDI Gridded Land Surface Metrics (LSM) L3 1km EASE-Grid, Version 2
## 5 GEDI Gridded Aboveground Biomass Density (AGBD) L4B 1km EASE-Grid, Version 2.1
## 6 MODIS/Terra+Aqua Land Cover Type (LC) Yearly L3 Global 500 m SIN Grid
## frequency resolution_meters
## 1 1 day 1000
## 2 Varies 70
## 3 Varies 70
## 4 One time 1000
## 5 One time 1000
## 6 1 year 500
# Check available bands for MODIS product MYD11A2
head(mt_bands("MYD11A2"))
## band description valid_range fill_value
## 1 Clear_sky_days Day clear-sky coverage 1 to 255 0
## 2 Clear_sky_nights Night clear-sky coverage 1 to 255 0
## 3 Day_view_angl View zenith angle of day observation 0 to 130 255
## 4 Day_view_time Local time of day observation 0 to 240 255
## 5 Emis_31 Band 31 emissivity 1 to 255 0
## 6 Emis_32 Band 32 emissivity 1 to 255 0
## units scale_factor add_offset
## 1 <NA> <NA> <NA>
## 2 <NA> <NA> <NA>
## 3 degree 1 -65
## 4 hrs 0.1 0
## 5 <NA> 0.002 0.49
## 6 <NA> 0.002 0.49
# Download MODIS surface temperature data
modis_temp <- mt_subset(
product = "MYD11A2",
lat = 46.0959,
lon = 15.1874,
band = "LST_Day_1km",
start = "2021-01-01",
end = "2022-12-31",
km_lr = 80,
km_ab = 80,
site_name = "Savinja-VelikoSirje",
internal = TRUE,
progress = FALSE
)
# Convert data to terra format
modis_temp_raster <- mt_to_terra(df = modis_temp)
# Plot surface temperature data (Kelvin)
plot(modis_temp_raster)
# Project the first layer (Jan 1, 2021)
selected_layer <- project(modis_temp_raster[[1]], crspost)
plot(selected_layer)
lines(catchment_area, col = "red")
# Extract temperature data for each date and convert to Celsius
temperature_data <- data.frame()
for (i in 1:nlyr(modis_temp_raster)) {
projected_layer <- project(modis_temp_raster[[i]], crspost)
projected_layer <- subst(projected_layer, 0, NA)
temp_kelvin <- extract(projected_layer, catchment_area, mean, rm)[2]
if (!is.na(temp_kelvin)) {
temp_celsius <- temp_kelvin - 273.15
new_row <- data.frame(Date = as.Date(names(modis_temp_raster)[i]), Temp_Celsius = temp_celsius, Temp_Kelvin = temp_kelvin)
colnames(new_row) <- c("Date", "Temp_Celsius", "Temp_Kelvin")
temperature_data <- rbind(temperature_data, new_row)
}
}
head(temperature_data)
## Date Temp_Celsius Temp_Kelvin
## 1 2021-02-10 0.7508639 273.9009
## 2 2021-02-18 13.5209944 286.6710
## 3 2021-02-26 13.3185523 286.4686
## 4 2021-03-06 7.8965254 281.0465
## 5 2021-03-22 16.3787452 289.5287
## 6 2021-03-30 21.1023461 294.2523
# Filter measured air temperature data
measured_temp <- main_data[main_data$Date >= as.Date("2021-01-01"), c("Date", "Airtemp")]
# Plot comparison between measured and MODIS temperature data
plot(
as.Date(measured_temp$Date), measured_temp$Airtemp,
type = "l", col = "blue", xlab = "Date", ylab = "Temperature (°C)",
main = "Temperature Comparison: Measured vs. MODIS",
ylim = c(-10, 35)
)
lines(as.Date(temperature_data$Date), temperature_data$Temp_Celsius, col = "red", type = "o", pch = 16)
legend(
"topright",
legend = c("Measured Data", "MODIS Data"),
col = c("blue", "red"),
lty = 1,
pch = c(NA, 16)
)
# Load ERA5 precipitation data (yearly, land domain)
# download_data("era5", getwd(), timestep = "yearly", domain = "land")
era5_precip <- rast("era5_tp_mm_land_195901_202112_025_yearly.nc")
# View basic properties and plot the first layer
era5_precip
## class : SpatRaster
## dimensions : 720, 1440, 63 (nrow, ncol, nlyr)
## resolution : 0.25, 0.25 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source : era5_tp_mm_land_195901_202112_025_yearly.nc
## varname : tp (Total precipitation)
## names : tp_1, tp_2, tp_3, tp_4, tp_5, tp_6, ...
## unit : mm, mm, mm, mm, mm, mm, ...
## time : 1959-01-01 to 2021-01-01 UTC
plot(era5_precip[[1]])
# Project catchment area to match ERA5 data
catchment_projected <- project(catchment_area, crs(era5_precip))
# Plot ERA5 data with catchment boundary overlay
plot(era5_precip[["tp_1"]], ext = c(14, 16, 45.5, 47))
lines(catchment_projected, col = "red")
# Extract grid cells within the catchment area
extracted_cells <- extract(era5_precip[[1]], catchment_projected, xy = TRUE, touches = TRUE, weights = TRUE)
head(extracted_cells)
## ID tp_1 x y weight
## 1 1 1236.111 14.625 46.375 0.35920
## 2 1 1234.035 14.875 46.375 0.71614
## 3 1 1278.850 15.125 46.375 0.75436
## 4 1 1233.905 15.375 46.375 0.43063
## 5 1 1129.915 15.625 46.375 0.00890
## 6 1 1233.344 14.875 46.125 0.16533
# Calculate weighted precipitation
weighted_precip <- sum(extracted_cells$tp_1 * extracted_cells$weight, na.rm = TRUE)
weighted_precip
## [1] 4197.993
annual_precip <- rep(NA, dim(era5_precip)[3])
# Compute average precipitation for each year
for (i in 1:dim(era5_precip)[3]) {
values <- extract(era5_precip[[i]], catchment_projected)[[2]]
annual_precip[i] <- mean(values, na.rm = TRUE)
}
# Plot annual ERA5 precipitation
time_series <- as.Date(time(era5_precip))
plot(
time_series, annual_precip, type = "o",
xlab = "Year", ylab = "Annual Precipitation (mm)",
main = "Annual ERA5 Precipitation for Savinja Catchment"
)
# Aggregate daily Areal Precipitation to annual totals
areal_precip_annual <- main_data %>%
group_by(Year) %>%
summarise(Areal_Precip = sum(Areal_Precip, na.rm = TRUE))
# Convert Year to Date format for plotting
areal_precip_annual$Year <- as.Date(paste0(areal_precip_annual$Year, "-01-01"))
# Filter ERA5 data from 2010 onward
time_vector <- as.Date(time(era5_precip))
start_index <- which(time_vector == as.Date("2010-01-01"))
time_vector_filtered <- time_vector[start_index:length(time_vector)]
annual_precip_filtered <- annual_precip[start_index:length(annual_precip)]
# Plot Areal and ERA5 Precipitation Comparison
plot(
areal_precip_annual$Year, areal_precip_annual$Areal_Precip, type = "o", col = "blue",
xlab = "Year", ylab = "Precipitation (mm)",
main = "Comparison of Areal and ERA5 Precipitation",
xlim = as.Date(c("2010-01-01", "2021-12-31")),
ylim = c(900, 1700),
xaxt = "n"
)
# Add ERA5 precipitation data
lines(time_vector_filtered, annual_precip_filtered, col = "red", lty = 1, type = "o")
axis(
1,
at = seq(as.Date("2010-01-01"), as.Date("2022-01-01"), by = "1 year"),
labels = format(seq(as.Date("2010-01-01"), as.Date("2022-01-01"), by = "1 year"), "%Y")
)
# Add legend
legend(
"topright", legend = c("Areal Precipitation", "ERA5 Precipitation"),
col = c("blue", "red"), lty = 1, pch = c(16, 16)
)